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A  User's  Guide  to  a  Computer  Program  for 
Harmonic  Analysis  of  Data  at  Tidal  Frequencies 

R.  E.  DENNIS  AND  E.  E.  LONG 
Oceanographic  Division,  Office  of  Marine  Surveys  and  Maps 


ABSTRACT.  This  report  describes  a  FORTRAN-IV  program  for  the  harmonic 
analysis  of  a  series  of  15  or  29  days  of  uniformly  spaced  tidal  data  (observations). 
The  program  is  usable  with  minor  modifications  on  any  computer  with  a  140K  mem- 
ory which  accepts  FORTRAN  input.  The  mathematical  basis  and  equations  for  the 
determina  tion  of  tidal  constants  from  observational  data  are  given.  This  report  should 
be  used  with  Manual  of  Harmonic  Analysis  and  Prediction  of  Tides,  C&GS  Special 
Publication  No.  98,  Revised  (1940)  Edition,  published  in  1941,  for  reference  to  for- 
mulas and  tables. 


It  was  recognized  in  ancient  time  that  tides  follow 
a  relatively  regular  cycle,  recurring  periodically. 
Since  any  periodic  motion  or  oscillation  can  be  re- 
solved into  components  consisting  of  simple  har- 
monic motions,  a  method  for  reducing  tidal  motions 
to  these  components  was  inevitable.  Lord  Kelvin 
devised  such  a  method  about  1867.  His  method 
followed  the  suggestions  and  speculations  made  by 
Laplace,  Young,  and  Airy  during  the  earlier  part  of 
the  19th  century,  but  full  credit  belongs  to  Kelvin 
for  making  such  analysis  practical.  In  the  late  19th 
century  and  early  20th  century,  William  Ferrel  and 
Rollin  A.  Harris  of  the  U.S.  Coast  and  Geodetic 
Survey  both  made  important  contributions  to  har- 
monic analysis  of  tides  and  interpretation  of  the 
results. 

In  1885,  Leland  P.  Shidy  of  the  U.S.  Coast  and 
Geodetic  Survey  designed  a  set  of  stencils  for  the 
"Standard  Method"  of  harmonic  analysis  of  obser- 
vational data  (Report  of  the  U.S.  Coast  and  Geo- 
detic Survey,  1893,  Vol.  1,  p.  108).  A  stencil  with 
apertures  was  fitted  over  tabulations  of  observed 
data  to  obtain  sums  for  computing  the  harmonic 
constants  of  the  tidal  constituents.  Only  after  this 
—  a  process  involving  computation  of  values  for 
about  20  coefficients  for  use  in  the  tide  equations — 
could  the  tide  prediction  begin.  Operational  har- 
monic analysis  was  continued  by  hand  until  the 


presently  used  computer  system  was  implemented 
in  1965. 

The  original  version  of  this  program  for  the 
standard  harmonic  analysis  of  tidal  data  was  com- 
pleted in  1965.  Since  then,  improvements  in  the 
method  and  refinements  in  programing  have  yielded 
a  program  that  can  analyze  the  observational  data 
for  a  station  in  as  little  as  1.5  seconds  on  a  CDC- 
6600  computer. 

BASIC  HARMONIC  EQUATION  AND 

DETERMINATION   OF  TIDAL   CONSTANTS 

FROM  OBSERVATIONAL  DATA 

For  analyzing  equally  spaced  short-period  data 
(15  days  or  29  days),  this  program  utilizes  the 
standard  Fourier  analysis  and  traditional  methods 
of  the  former  Coast  and  Geodetic  Survey  described 
by  Schureman  (1941)  with  extracts  from  Doodson 
(1924)  and  Harris  (1897)  reports  on  tidal  analysis. 

Harmonic  analysis  of  tidal  data  entails  three 
basic  processes: 

1.  Initial  separation  of  the  tidal  constituents  from 
the  data 

2.  Orientation  of  the  constituent  tides  with  the 
astronomical  elements 

3.  Elimination  from  each  constituent  sought  of 
effects  of  other  tidal  constituents 


The  equation 

h  =  Ho  +^2fn  Hn  cos  [ant  —  (*„  —  [V0  +  u]n)] 
n 

(1) 
describes  the  height  of  tide  at  any  time  for  which 
Ho — mean  value  of  tide  for  observations 

1    N~l 

Ho  —  -z}  2-i  hi 

N    i=0 

t — time  reckoned  from  some  arbitrary 
point 

/„ — node  factor  of  constituent 

Hn — mean  amplitude  of  the  constituent 

an — speed  of  constituent 

Kn — epoch  of  the  constituent  at  t  =  0 
for  period  of  observation 

(Fo+u)„ — value  of  argument  of  constituent  at 
t  =  0  for  period  of  observation 

n — any  particular  constituent  being 
computed 

iV — total  number  of  data  values  used 
in  analysis 

The  factor  Hn  is  the  mean  amplitude  of  constituent, 
n,  for  an  entire  nodal  period,  whereas  the  adjust- 
ment using  the  node  factor /„,  denoted 


.,  (N  —  1)5  where  N8  =  2w  provided  that 


J^-n  Jn  **-n 


(2) 


is  the  amplitude  pertaining  to  a  particular  time. 
The  amplitudes  Rn  are  derived  from  observational 
data. 

Although  the  constant  Kn  is  used  ultimately  to 
characterize  the  constituent,  the  quantity  k„  — 
{V0-\-u)n  may  be  designated 


$ n    =    Kn  —    (V0-\-u)n 


(3) 


where  f„  is  derived  from  the  data.  £"„  represents  the 
phase  of  the  constituent  at  t  =  0  for  a  particular 
series  of  observations. 

To  use  the  tidal  equation  in  harmonic  form  we 
follow  Schureman's  (1941)  discussion  of  the  Fourier 
series  technique.  The  curve  described  by 

h  =  H0  +  Ci  cos  0+  d  cos  20  +  .  .  .  Ck  cos  &0 
+  Si  sin  0+  S2  sin  20  +  .  .  .  Sz  sin  Z0 

(4) 

will  pass  through  the  N  coordinates  given  by  ordi- 
nates  ho,  hi,  hi,  hz,  .  .  .,  Ajv-i  and  abscissas  0,  5,  25, 


2     y^1 
Cp  =  -r7  2-i    hi  cos  ipd 

J\    ;=o 


and 


N—  1 


sD  = 


Tf   /Li    hi  sin  iph 


(5) 


(6) 


Since  the  limits  of  k  and  I  are 

N 


for  N  even 


lim  k  = 
lim  / 


-  1 

AT-1 


for  TV  odd,  lim  k  =  lim  I  = 

z 

N 
the  limit  of  p  is  correspondingly  always  less  than  — . 

When  0  represents  the  constituent  day,  Ci  and  Si 
corresponding  to  p  =  1  represent  the  diurnal  con- 
stituents, and  Ci,  Si  corresponding  to  p  =  2  repre- 
sent the  semidiurnal  constituents.  For  tidal  work 
values  of  p  greater  than  8  are  seldom  used,  so  the 
limit  of  N/2  has  no  significance. 

Using  the  notation 

a  =  speed  of  constituent  sought  in  degrees 
per  solar  hour 

5  =  length  of  sampling  period  in  solar  hours 

the  coefficient  equations  may  be  written  more  ex- 
actly for  real  data: 


C*  =  !f§  hiCos[(i-l)[^a8p)\ 


(7) 


(8) 


From  the  Fourier  coefficients  in  equations  (7)  and 
(8)  desired  harmonic  constants  are  defined: 


r^tan-d) 

R'(A)  =  (SP2  +  Cp2)1'2 


(9) 
(10) 


The  quantity  $'{A)  is  the  phase  at  t  =  0,  and  R'{A) 
the  amplitude  for  the  constituent  during  the  par- 
ticular series  of  observations.  Both  symbols  are 
primed,  designating  that  the  effects  of  other  con- 
stituents have  not  been  removed.  For  each  constit- 
uent the  quantity  p  in  the  coefficient  equations 
equals  the  coefficient  of  the  hour  angle  of  the  mean 
sun,  T,  in  the  equilibrium  argument  of  that  con- 
stituent. 


Subroutine  FORAN  in  the  program  does  the  entire 
Fourier  summation  for  each  desired  constituent.  R' 
and  f'  are  determined  from  the  original  series  for 
five  major  constituents  M2,  S2,  N2,  Oi,  and  Ki  and 
the  harmonics  M4,  M6,  Ms,  S4,  and  S6.  In  the  case 
of  a  15-day  series  of  observations,  the  N2  constit- 
uent is  approximated  from  the  constants  for  M2. 

Astronomical  Adjustments 

Modifications  of  R'  and  f '  must  be  made  to  reduce 
these  constants  from  the  particular  time  of  obser- 
vation to  mean  values  for  an  entire  nodal  period. 
The  epochs  of  each  constituent  (F0  +  u)n  must  be 
computed  to  adjust  f '„  to  k  n  and  node  factors  fn 
determined  to  resolve  R'n  to  H'n  for  the  entire 
nodal  period. 

Calculation  of  Orbital  Functions 

Prior  to  either  adjustment,  certain  orbital  elements 
and  functions  of  these  elements  must  be  calculated. 
These  are: 

Primary  Orbital  Elements 

s — mean  longitude  of  the  moon 
p — mean  longitude  of  lunar  perigee 
h — mean  longitude  of  the  sun 
pi — mean  longitude  of  solar  perigee 
N — longitude  of  ascending  lunar  node 
i — inclination  of  moon's  orbit  to  plane  of 

ecliptic 
co — obliquity  of  the  ecliptic 

Secondary  Orbital  Functions 

I — obliquity  of  the  moon's  orbit 

v — longitude  in  the  celestial  equator  (right  as- 
cension) of  the  moon's  intersection  with  the 
equator 

£ — longitude  in  the  moon's  orbit  of  the  moon's 
intersection  with  the  celestial  equator 

v' — a  function  of  the  moon's  orbit  (this  term 
appears  in  the  argument  of  lunisolar  con- 
stituent Ki) 
1v"  —  another  function  of  the  moon's  orbit  (this 
term  appears  in  the  argument  of  lunisolar 
constituent  K2) 

P — mean   longitude   of  lunar   perigee   reckoned 
from  lunar  intersection  in  the  celestial  equator 

Each  of  the  orbital  elements  must  be  calculated 
for  the  beginning  of  the  series  before  the  epochs  can 
be  computed.  The  subroutines  ASTRO  and  ORBIT 
compute  these  elements.  Formulas  for  their  determi- 


nation used  in  the  program  and  found  in  table  1 
are  compiled  from  Newcomb  (1912)  and  Doodson 
(1924).  Detailed  discussion  of  these  elements  can  be 
found  in  these  references  and  in  other  works  on 
celestial  mechanics. 

Computation  of  the  secondary  orvital  functions 
are  based  on  the  assumption  i  and  co  are  constant. 
This  assumption  may  be  reliable  for  as  much  as  4 
centuries  prior  to  or  after  the  year  1900,  as  co  is  de- 
fined by  the  function 

co  =  23.45229444°  —  0.0130125°T 

—  0.000001638°r2 

and  varies  only  0.1  degree  in  4  centuries.  (T  defined 
in  table  1.)  Formulas  for  the  secondary  orbital  func- 
tions are  given  in  table  2,  based  on  Baird  (1886) 
and  Schureman  (1924). 

Computation  of  the  Equilibrium  Arguments 

Equilibrium  arguments  ( V  +  u)  are  combinations 
of  the  lunar  and  solar  elements  in  the  tidal  equations 
derived  from  lunar  and  solar  motion  and  are  de- 
noted by  the  symbol  E(A).  They  serve  to  identify 
the  constituent  and  determine  its  speed  and  period, 
fixing  times  of  extrema  for  the  corresponding  tidal 
force.  When  referring  the  argument  E  to  the  time 
at  the  beginning  of  a  particular  series,  the  argument 
becomes  Vq  +  u.  The  basic  equations  for  computing 
the  equilibrium  arguments  are  taken  from  Schure- 
man (1924).  Equilibrium  arguments  for  (MK)j, 
(2MK)3,  (MN)4,  (MS)4,  (2SM)2,  Mf,  MSf,  Mm,  Sa, 
and  Ssa  are  omitted  because  they  have  negligible 
effect  upon  the  maxima  and  minima  for  a  short - 
period  analysis.  M3  and  Si  are  also  omitted  because 
their  effects  cannot  be  determined  from  a  short 
period  analysis. 

Computation  of  Node  Factors 

The  ratio  of  the  true  obliquity  of  the  moon's  orbit 
for  any  particular  value  of  J  to  the  mean  value  of 
the  obliquity  over  an  entire  nodal  period  is  called 
the  node  factor.  This  factor  applied  to  the  mean 
constituent  amplitude  for  a  nodal  period  yields  the 
amplitude  for  the  constituent  at  that  particular 
value  of  /.  Conversely,  the  reciprocal  of  the  node 
factor  applied  to  the  tidal  coefficients  determined 
from  a  particular  series  reduces  them  to  a  mean 
value  for  the  constituent  for  all  positions  of  the 
lunar  node.  Equations  used  for  computation  of  the 
node  factors  were  extracted  from  Harris  (1897)  and 
Schureman  (1941)  and  are  shown  in  table  3. 


Table  1. — Formulas  for  computation  of  orbital  elements 
Formulas  for  adjustment  of  orbital  elements  to  beginning  of  century: 


P 
h 


=  270.43742222°  +  307.892°T  +  0.002525°T2  +  0.00000189°T3  +  s0z 

=  334.32801944°  +  109.0322055°T  -  0.01034444°T2  —  0.0000125°T3  +  p0. 

=  279.69667778°  +  0.768925°T  +  0.0003025°T2  +  h0z 

=  281.22083333°  +  1.719175°T  +  0.00045278°T2  +  0.00000333°T3  +  ploz 


N'  =  259.18253333°  -  134.1423972°T  +  0.002 10556° T2  +  0.00000222° +  T3  +  NoZ 

Formulas  for  adjustment  of  orbital  elements  to  beginning  of  observations: 

[For  beginning  of  the  series] 

s  =  s'  +  129.38482032  [Y-C]  +  13.176396768  [Ds  +  X]  +  0.549016532    [grbs] 

p  =  p>  +  40.66246584    [Y-C]  +  0.111404016    [Ds  +  X]  +  0.004641834    [grbs] 

h  =  h'  -  0.238724988    [Y-C]  +  0.9856473288  [Ds  +  X]  +  0.0410686387  [grbs] 

pi  =  pi  +  0.1717836        [Y-C]  +  0.000047064    [Ds+X]  +  0.000001961    [grbs] 

[For  middle  of  series] 

p  =  p>   +  40.66246584    [Y-C]  +  0.111404016    [Dm  +  X]  +  0.004641834    [grms] 
N  =  N'  +  19.328185764  [Y-C]  -  0.0529539336  [Dm+X]  -  0.0022064139  [grms] 

s0,  Po,  h0,  pu,  N0 — Speeds  of  the  elements  (degrees/solar  day) 

s',  p' ,  h',  pi,  N' — Values  of  the  elements  at  the  beginning  of  any  century 

s,  p,  h,  pi,  N — Values  of  the  elements  at  any  given  time 

T — Nearest  whole  number  of  Julian  centuries  (36525  days,)  from  Greenwich  Mean 
Noon,  December  31,  1899,  on  the  Gregorian  Calendar  to  the  beginning  of  the 
century  for  observations 

Y — The  exact  year  in  which  observations  are  recorded 

C — The  century  during  which  observations  are  recorded 

z — A  factor  representing  the  time  difference  (in  solar  days)  formed  between  the 
Gregorian  calendar  date  and  the  time  T  as  calculated  in  Julian  centuries 

Da — The  day  of  the  year  during  which  observations  begin 

Dm — The  day  of  the  year  upon  which  the  midpoint  of  the  series  falls 

X — Correction  for  leap  years 

X  =  0.25  [F— (c+1)] 

truncated  to  no  decimals 

grbs — Greenwich  hour  at  the  beginning  of  the  series 

grms — Greenwich  hour  at  the  middle  of  the  series 


Table  2. — Formulas  for  secondary  orbital  functions 

I  =  cos-1  [0.9136949  -  0.0356926  cos  N]* 

v  =  sin-1  [0.0897056  sin  iV/sin  I]* 

_j  I"  (0.206727  sin  N)  (1  -  0.0194926  cos  N)         >* 

*  ~  tan     Lo.9979852  +  0.206727  cos  N  —  0.0020148  cos  2NJ 

.  T  •  I  ,    0.334766^1** 

"  = tan  Lsin ' +  Vcos " +  ~^2r) J 

i  T  ■    o         (       o     ,   0.0726184M** 
2j>"  =  tan-1    sin  2^  -J-  I  cos  2^  + .  ) 

P  =  [p  (for  middle  of  series)  —  £]* 
u  of  Z,2  =  2  ({  -  v)  -  R  ±  180°* 
where  R  =  tan"1  [sin  2P  /  { (cos2  0.5/  /  6  sin2  0.57)  —  cos  2P}  ] 

uofMi  =  U-v)  +(?  +  90°* 
where  Q  =  tan"1  [  { (5  cos  I  —  1)  /  (7  cos  I  +  1) }  tan  P] 

*Schureman  (1924) 
**Baird  (1886) 

Table  3. — Formulas  for  calculation  of  node  factor  reciprocals 

F(A)-/(A)  =  1.000 

/(A)  =  node  factor  of  constituent  A. 

F(M2)  =  (0.91544)  ^-  (cos4  0.5  I") 

F(S2)  =  F(R2)  =  F(T2)  =  F(P0  =  F(S4)  =  F(S«)  =  1.000 

F(Oi)  =  (0.37988)  -f-  (sin  7  cos2  0.5  I) 

F(Ki)  =  (0.8965  sin2  2/  +  0.6001  sin  21  cos  k  +  0.1006)"1/2 

F(K2)  =  (19.0444  sin4  7  +  2.7702  sin2  7  cos  2v  +  0.0981)"1/2 

F(L2)  =  F(M2)  X  [1  -  (12  sin2  0.5  I  cos  2P)  +  (cos2  0.5  I)  + 

(36  sin4  0.5  I)  -i-  (cos4  0.5  J)]~1/2 

F(J0  =  (0.72137)  -f-  sin  21 
F(Mi)  =  F(Oi)  X  (2.310  +  1.435  cos  2P)~1/2 
F(OO)  =  (0.016358)  -=-  (sin  J  sin2  0.5  I) 

F(N2)  =  F(2N)  =  F(X2)  -  Ffx2)  =  F(v2)  =  F(M2) 

F(Q0  =  F(2Q)  =  F(P0  =  F(O0 
F(M4)  =  F2(M2) 
F(M«)  =  F3(M2) 
F(Mg)  =  F4(M2) 


Inference  of  R',  k',  f'  for  Secondary 
Constituents 

Following  the  discussion  of  Schureman,  the  re- 
lations existing  between  observational  constants  of 
similar  constituents  are  assumed  the  same  as  those 
existing  between  the  corresponding  theoretical 
values.  Since  values  for  only  10  constituents  are 
sought  (M2,  S2,  N2,  Oi,  Ki,  M4,  M6,  M8,  S4,  S6),  the 
other  constituents,  although  minor  in  effect,  must 
be  approximate  before  the  process  of  elimination 
can  be  applied  to  the  major  constituents.  In  the 
program  this  is  done  using  the  formulas  given  by 
Schureman  (1941). 

Elimination  of  Perturbations  From 
Secondary  Constituents 

The  disturbing  effects  upon  S2  by  K2  and  T2  and 
on  Kx  by  Pi  may  be  considerable  in  a  short  series 
owing  to  the  small  differences  in  frequencies  of  the 
respective  constituents.  To  correct  for  this,  account 
is  taken  of  the  phase  displacement  and  augmenta- 
tion of  amplitude  in  the  major  constituents  S2  and 


Ki  by  the  respective  minor  constituent  disturb- 
ances. This  process  is  described  by  Schureman 
(1941)  and  performed  by  the  computer  using  the 
imbedded  data  from  tables  21-26  from  the  Manual 
of  Harmonic  Analysis  and  Prediction  of  Tides. 

The  final  process  of  elimination  frees  the  constit- 
uents from  the  disturbing  effects  of  the  other  con- 
stituents. The  constituents  sought  are  M2,  N2,  S2, 
Oi,  Ki.  Only  constituents  with  the  same  coefficient 
for  the  factor  T  in  their  equilibrium  argument  are 
considered  interfering,  that  is,  effects  of  diurnal 
upon  semidiurnal  and  vice  versa  are  negligible. 

Complete  discussion  of  constituent  interference 
is  found  in  Schureman  (1941).  The  program  com- 
pletes the  elimination  process,  using  tables  prepared 
from  equations  389,  390  from  Schureman  (1941, 
table  29). 

Harmonic  Constants 

Inference  of  secondary  constituents  from  primary 
constituents  utilizes  the  ratios  derived  by  Schure- 
man (1941)  and  shown  in  table  4. 


Table  4. — Formulas  for  inference  of  secondary  constituents 

Diurnal  constituents 

//(JO     =  0.079  tf(Oi);  *(Ji)     =  k(Kx)  +  0.496  [k(Kx)  -  *(Oi)] 

#(Mi)  =  0.071  jff(Oi);  «(Mi)  =  k(Ki)  -  0.500  [K(Ki)  -  k(Oi)] 

#(00)  =  0.043  i/(Oi);  k(OO)  =  k(Ki)  +  1.000  [*(Ki)  -  *(Oi)] 

f/(Px)    =  0.331  //(KO;  k(Pi)    =  k(Kx)  -  0.075  [k(Ki)  -  k(Oi)] 

#(Qi)    =  0.194  #(00;  k(Qi)    =  k(Ki)  -  1.496  [«(Ki)  -  k(Oi)] 

H(2Q)   =  0.026  i/(Oi);  k(2Q)    =  k(Kx)  -  1.992  [k(K0  -  K(Ox)] 

H(P1)     -  0.038  //(Ox);  «(Pl)     =  K(KX)  -  1.429  [k(Ki)  -  k(Oi)] 

Semidiurnal  constituents 

ff(K2)    =  0.272  H(Sf);    k(K2)  =  K(S2)    +  0.081  [/c(S2)    -  /c(M2)] 

H(U)    =  0.028  tf(M2);  «(L2)  =  k(S2)    -  0.464  [k(S2)    -  /c(M2)] 

=  0.143  H(Ni);  =  «(M2)  +  1.000  [k(M2)  -  /c(N2)  ] 

i/(N2)    =  0.194  H(M2);  «(N2)  =  «(S2)    -  1.536  [«(S2)    -  X(M2)] 

H(2N)  =  0.026  H(M2);  *(2N)  =  *c(S2)    -  2.072  [k(S2)    -  /c(M2)] 

=  0.133  H(N2);  =  k(M2)  -  2.000  [*(M2)  -  k(N2)  ] 

i/(R2)    =  0.008  ff(S2);    /c(R2)  =  k(S2)    +  0.040  [/c(S2)    -  k(M2)} 

H(T2)    =  0.059  H(S2);     k(T2)  =  k(S2)    -  0.040  [k(S2)    -  «(M2)] 

J7(A2)     =  0.007  i/(M2);  k(X2)  =  k(S2)    -  0.536  [k(S2)    -  /c(M2)] 

H(n2)     =  0.024  H(M2);  *0*,)  =  k(S«)    -  2.000  [«(S2)    -  k(M2)] 

J7(y2)     =  0.038  i/(M2);  k(v2)  =  k(S2)    -  1.464  [k(S2)    -  k(M2)] 

=  0.194  H(N2);  =  k(M2)  -  0.866  [«(M2)  -  K(N2)  ] 


USERS'  GUIDE  TO  PROGRAM 

The  present  program  is  written  in  FORTRAN- 
IV,  executable  with  minor  adjustments  on  any 
compatible  machine  having  a  140K  memory  and 
access  to  arcsine  and  arccosine  systems  functions. 
Computing  time  is  approximately  1.5  seconds  on 
the  CDC -6600. 

Either  a  vector  (polar  form)  or  scalar  variable 
may  be  analyzed.  For  vector  series,  the  program 
allows  either  a  major-minor  axis  analysis  or  a  north- 
east component  approach.  No  data  series  may  ex- 
ceed 7,000  terms  without  redimensioning  in  the 
program,  and  no  series  of  other  than  15  or  29  days 
of  uniformly  spaced  data  can  be  analyzed. 

The  program  accepts  input  via  magnetic  tape  or 
punched  cards  in  any  format  with  the  restriction 
that,  for  vectors  with  magnitude  and  direction  in 
the  same  record,  the  angles  must  precede  the  am- 
plitudes in  the  record.  For  vectors  specified  by  one 


file  of  amplitudes   and   one   file  of  directions,   the 
amplitude  file  must  be  read  first. 

Output  comprises  the  mean  amplitudes  and  phases 
of  26  tidal  constituents. 

Explanation  of  Control  Cards 

The  order  of  the  control  cards  (fig.  1)  input  to  the 
machine  is  of  utmost  importance.  The  discussion  of 
each  card  or  set  of  cards  is  in  the  order  in  which 
they  are  read.  The  master  job  control  card  appears 
only  once,  whereas  all  other  cards  appear  for  each 
new  data  series. 

Master  Job  Control  Card 

Immediately  following  the  program  deck  must 
be  a  card  containing  the  quantity  NJ.  This  specifies 
the  number  of  different  jobs  to  be  done. 

READ  221,  NJ 
221  FORMAT  (15) 
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Figure  1. — Order  of  control  cards.  If  data  are  on  cards,  data  deck(s)  follow  the  format 
specification  card  and  the  correction  cards  are  not  applicable. 


Identification  Cards 

These  two  cards  follow  the  master  job  control 
card  and  are  read  by  the  statement. 

READ  808,  IDENS 
808  FORMAT  (9A8/9A8) 

The  purpose  of  these  cards  is  to  introduce  any  de- 
sired alphanumeric  information  to  be  used  as  iden- 
tifications. Contents  of  these  two  cards  will  be 
printed  verbatim  immediately  preceding  final  results. 
Only  columns  1-72  of  each  card  will  be  printed. 

Program  Control  Card 

The  program  card  specifies  the  operations  to  be 
performed  on  the  data  series,  the  mode  of  input  for 
data  series,  and  any  adjustments  which  are  to  be 
made  to  the  data.  The  variables  are  read  by  the 
statement 

READ  76,  AZI,  LSORT,  JJ,  N,  MAP,  K, 
CVAR,  INC,  NSPH,  TFAC,  EAZI 
76  FORMAT  (F5.0,  515,  F5.0,  215,  F10.8,  F5.0) 

The  variables  listed  have  the  following  significance: 

AZI — Approximate  mean  flood  direction  (F5.0) 

LSORT — Specifies  when  a  sort  of  the  input  data 
is  required  (15) 

LSORT  =  0.  No  sort 

LSORT  =  1.  Sort  to  be  performed 

JJ — Indicates  the  type  of  sort  desired  (15)   (See 
subroutine  SORT) 

N — The  number  of  data  points  in  the  series  (15) 

MAP — Controls  input  mechanism  and  format. 
MAP  is  determined  partially  by  whether 
the  data  are  vector  or  scalar.  (15) 

MAP  =  1.  Scalar  quantity,  card  input  with 
format  AQ.  (See  format  specification 
card) 

MAP  =  2.  Magnetic  tape  input  for  vector 
quantity  with  format  AQ.  Angles  must 
precede  magnitudes  in  the  record 

MAP  =  3.  Card  input  of  vector  data  with 
magnitudes  of  format  AQ  and  direc- 
tions of  format  RQ. 

MAP  =  4.  Tape  input  of  vector  data  accord- 
ing to  format  AQ  and  tilt  corrections 
from  tape  (if  available)  according  to 
format  BQ.  Use  MAP  =  4  only  for 
making  corrections  to  data  series.  (See 
subroutine  EDAT) 


K — Indicates  whether  analysis  by  components  is 
desired  (15) 

K  =  1.  Major  axis  analysis  only  (use  K  =  1  for 
scalar) 

K  =  2.  Minor  axis  analysis  only 

K  =  3.  Major  and  minor  axis  analysis 

K  =  4.  North  and  east  component  analysis.  In 
this  case  AZI  =  0.00 

CVAR — Compass  variation  of  the  vicinity.  CVAR 
is  negative  when  West  and  positive  when 
EAST.  (F5.0) 

INC — Indicate  if  tilt  corrections  are  to  be  applied. 
(See  subroutine  TILT)  (15) 

INC  =  0.  No  tilt  corrections 

INC  =  1.  Tilt  corrections  to  be  applied 

NSPH  —  Number  of  samples  of  data  per  mean 
solar  hour  (15) 

TFAC — Time  correction  factor  for  a  slow  or  fast 
chronometer  in  meters.  (F10.8)  For  no 
time  correction  TFAC  =  1.000000.  To 
determine  TFAC  use  the  expression 

—,„  .  „  _    True  time  covered  by  observations 
Time  reckoned  by  the  meter  clock 

CAUTION:  If  NSPH  or  TFAC  is  zero,  the  program 
will  terminate  immediately. 

EAZI  —  approximate  ebb  azimuth  (F5.0) 

Astronomical  Parameter  Card 

The  card  specified  in  the  statement 

READ  2038,  XYER,  MONTH,  DAY,  STT, 
TM,  GONL 
2038  FORMAT  (F5.0,  15,  F5.0,  F5.2,  2F6.2) 

inputs  information  which  is  pertinent  to  the  partic- 
ular times  of  the  series  of  observations. 

XYER — Year  in   which   the   series   is   observed 

(F5.0) 

MONTH— Month  during  which  series  starts  (15) 
DAY — Day  on  which  series  starts  (F5.0) 

STT — Time  in  decimal  hours  at  which  series  starts 

(F5.2) 

TM— Time  meridian  to  which  station  is  referred 

(F6.2) 

GONL— Longitude  of  the  station  (F6.2) 


CAUTION:  The  convention  of  west  longitude  being 
positive  and  east  longitude  being  negative  is  fol. 
lowed  in  entering  the  parameters  TM  and  GONL- 


Format  Specification  Card 

This  card  allows  variability  in  the  format  of  data 
on  magnetic  tape  or  cards.  This  card  is  read  by  the 
statement 

READ  2039,  AQ,  BQ 
2039  FORMAT  (4A8,  3A8) 

The  maximum  field  length  of  AQ  is  32  characters 
and  of  BQ  is  24  characters.  The  format  specification 
BQ  must  always  begin  in  column  33  of  this  card. 
Parenthesis  must  be  punched  and  counted  as  part 
of  the  field  length. 

Limitations  are  according  to  the  variable  MAP: 

MAP  =  1.  Assign  AQ  a  card  format;  do  not 
assign  BQ. 

MAP  =  2.  Assign  AQ  a  tape  format;  do  not 
assign  BQ.  Angle  of  vector  must  precede 
the  magnitude  in  the  record. 

MAP  =  3.  Assign  AQ  and  BQ  formats  for  card 
reading.  AQ  must  be  magnitudes  of  the 
vectors  and  BQ  must  be  directions. 

MAP  =  4.  AQ  and  BQ  must  be  as  follows: 
AQ  =  (60X,  F3.0,  24X,  F4.2) 
BQ  =  (110X,  F3.0) 

This  is  necessary  for  tilt  corrections. 


Correction  Control  Card 

If  data  are  on  cards,  data  decks  immediately  fol- 
low the  format  specification  card.  For  data  read 
from  magnetic  tape,  corrections  can  be  made  using 
the  correction  control  card. 

READ  222,  MX,  LRD,  LX 
222  FORMAT  (315) 

MX — Total  number  of  records  which  will  be  read 
to  attain  desired  series 
MX  =  N  +  LRD 

LRD — Number  of  records  at  beginning  of  tape 
to  be  skipped 

LX — Number  of  corrections  to  be  read  from  cards 
and  substituted  in  the  series  for  those  values 
found  on  the  magnetic  tape 


Correction  cards  follow  the  correction  control  card 
and  contain  the  interval  number  of  the  point  to  be 
changed  and  the  desired  angle  and  magnitude  of 
the  vector  at  that  point.  The  statement 

READ  555,  (ICOR  (LY),  CD  (LY), 
CV  (LY),  LY  =  1,  LX) 
555  FORMAT  (15,  F10.0,  F10.2) 

reads  these  parameters. 

Explanation  of  Subroutines 

The  subroutines  called  in  program  CURAN  can  be 
grouped  into  three  categories.  The  mathematical  anal- 
ysis and  astronomical  adjustments  are  performed  in 
four  basic  subroutines— FOR  AN,  ORBIT,  ASTRO, 
SASTR.  The  subroutine  DAYXX  is  a  utility  routine 
used  for  adjusting  the  time  of  beginning  of  series  to 
astronomical  reference  frame  (Greenwich  Time  Me- 
ridian). Operations  peculiar  to  tidal  currents  are 
performed  by  subroutines  AZIM,  DETAZ,  and 
TILT. 

Utility  routines  which  are  called  or  in  some  cases 
are  optional  include  SORT,  EDAT,  MEAN,  TERPO, 
TEMP,  FITAN,  ONEPI,  and  TWOPI. 


Subroutine  FOR  AN  (SER,  SP,  RES,  ZEP,  N,  M) 

This  subroutine  is  designed  to  determine  the 
Fourier  coefficients  of  the  data  (SER)  for  the  first 
approximations  of  phases  (ZEP)  and  amplitudes 
(RES)  of  the  various  constituents.  The  angular  dis- 
placement per  sampling  period  (SP)  of  the  required 
constituent  is  computed  in  the  subroutine  according 
to  the  expression: 


SP  = 


(Speed  of  the  constituent  in  degrees  per 

mean  solar  hour) 

(Number  of  samples  per  hour  X  180) 


The  numbers  N  and  M  represent  the  number  of  data 
points  and  the  number  of  harmonics  sought,  respec- 
tively. The  subroutine  uses  equations  (289)  and  (292) 
from  Schureman  to  obtain  the  Fourier  coefficients. 


Subroutine  SASTR 

This  subroutine  calls  as  input  the  reciprocals  (CF) 
of  the  node  factors  of  the  various  constituents  com- 
puted in  the  main  program  according  to  equations 

(73)-(80),  (141)-(150),  (207),  (215),  and  (235).  The 
input  also  comprises  the  node  factor  for  K2  (CKKF), 
the  mean  longitude  of  the  sun  (SL),  the  arguments 


for  v',  v"  (VP,  VPP),  and  the  mean  longitude  of  solar 
perigee  (PS).  Tables  21-26  from  Schureman  are  em- 
bedded in  the  program.  The  following  effects  of 
secondary  constituents  are  computed: 

ACS — Acceleration  in  S2  from  K2  and  T2 

APS — Augmentation  of  amplitude  of  S2  from  K2 

and  T2 

ACK — Acceleration  in  Ki  from  Pi 

APK — Augmentation  of  amplitude  of  Ki  from  Pi 

This  subroutine  calls  subroutine  TERPO  to  inter- 
polate the  appropriate  tables. 


Subroutine  ASTRO  (XYER,  DAYB,  DAYM, 
GRBS,  GRMS,  JOBX) 

COMMON/COSTX/CXX(30),  OEX(5) 

This  subroutine  uses  input  from  subroutine 
DAYXX  (DAYB,  DAYM,  GRBS,  GRMS)  with  the 
input  value  XYER,  the  year  in  which  observations 
begin,  to  compute  the  orbital  elements  s,  />,  h,  pu 
N,  I,  v,  E,  /,  and  2v  "  (table  1)  relative  to  the  period 
of  observation.  The  subprogram  calls  subroutine 
ORBIT  to  compute  all  values  for  the  beginning  of 
the  century,  then  adjusts  each  to  the  time  of  the 
observed  series.  The  array  CXX  are  the  appropriate 
values  of  the  orbital  elements.  JOBX  is  the  number 
of  days  in  the  data  series. 


Subroutine  ORBIT  (XCEN,  XSX,  SPX,  SHX, 

XP1X,  XNX,  OEX,  T,  XYER,  NNN) 

This  subroutine  uses  the  value  XYER  to  compute 
orbital  elements  (table  1)  for  the  beginning  of  the 
century  during  which  observations  were  taken.  The 
elements  are  computed  according  to  the  expressions 
given  by  Newcomb  and  Schureman  displayed  in 
table  1.  Values  are  computed  in  array  OEX  and 
transferred  to  the  variables  XSX,  XPX,  XHX, 
XP1X,  XNX  for  convenience.  XCEN,  T,  NNN 
are  all  internal  variables. 


Subroutine  TERPO 

This  subprogram  is  a  utility  program  designed  to 
interpolate  a  table  of  figures  in  two  directions. 


new  array  consists  of  the  first  value  and  every  JJth 
value  thereafter  of  the  original  series.  KK  is  internal 
and  dimensions  the  new  array  VX. 

Subroutine  EDAT  (MX,  VIN,  DIN,  MA,  LX, 

ICOR,  CD,  CV,  LRD,  XTEMP,  INC) 

This  routine  is  available  only  to  MAP  =  4.  It  is 
designed  for  cases  where  polar  vectors  are  input  as 
data  from  magnetic  tape,  but  either  simple  substi- 
tutions or  tilt  corrections  need  to  be  applied  to  the 
data  series. 

MX-MA  +  LRD 

VIN — Magnitudes  (corrected) 

DIN — Directions  (corrected) 

MA — Total  number  of  values  desired  in  final 
data  series 

LX — Number  of  substitutions  and  inferences 
to  be  made  within  the  series  or  at  end 
of  series 

ICOR — Interval  number  in  the  list  MX  which 
is  to  be  corrected 

CD — Angle  to  be  substituted  for  the  value 
actually  on  tape 

CV — Magnitude  to  be  substituted  for  the 
value  actually  on  tape 

LRD — Number  of  records  on  tape  to  be  skip- 
ped 

XTEMP — Temporary  storage  for  TILT  values 

INC — Integer  indicating  when  tilt  corrections 
are  to  be  applied 

CAUTION:  If  a  magnitude  or  direction  for  a  given 
interval  need  correcting,  both  magnitude  and  direc- 
tion and  interval  number  must  be  punched  on  a 
card.  The  program  substitutes  the  entire  record.  For- 
mat for  correction  cards  bearing  ICOR,  CD,  CV  is 
(15,  F10.0,  F10.2).  The  number  of  substitutions  used 
may  not  exceed  300. 

If  an  inference  of  data  is  required  to  extend  the 
series  to  the  required  length  of  15  or  29  days,  the 
magnitudes  and  directions  must  be  punched  on  cards 
with  the  formats  (24F3.2)  and  (24F3.0),  respectively. 
These  cards  then  follow  any  substitution  cards  in 
the  data  deck. 


Subroutine  SORT  (SORX,  VK,  KK,  JJ) 

This   is   designed   to   sort   values   from   an   array 
(SORX)  forming  a  new  array  (VX)  such  that  the 


Subroutine  XMEAN 

This  is  a  utility  routine  to  compute  the  arithmetic 
mean  of  an  array. 
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Subroutine  DAYXX 


Subroutine  DETAZ 


Computes  the  day  of  the  year  and  the  Greenwich 
hour  for  the  beginning  and  middle  of  the  data  series. 

Subroutine  TEMP 

Provides  temporary  storage  for  any  array  so  that 
operations  can  be  performed  on  the  array  without 
loss  of  the  original  series. 

Subroutine  FITAN 

This  subprogram  finds  the  angle  for  a  given  value 
of  the  arctangent  function. 


This  routine  is  designed  to  compute  a  mean  major 
axis  azimuth  for  a  given  approximate  azimuth  which 
is  input  on  the  master  job  contol  card. 

Subroutine  TILT 

This  subroutine  is  designed  for  correcting  errone- 
ous tidal  current  velocities  for  instrumental  tilt.  It 
should  be  applied  only  to  instruments  having  Sa- 
vonius  rotor  current  sensors  with  tilt  meters.  If  the 
tilt  reading  exceeds  35  degrees,  a  manual  inference 
will  be  required. 


Subroutine  ONEPI 

This  subroutine  is  designed  to  prevent  the  differ- 
ence between  two  angles  from  exceeding  180  degrees 
or  becoming  negative. 

Subroutine  TWOPI 

This  subroutine  places  a  given  angle  in  the  proper 
quadrant. 

Subroutine  AZIM  (DX,  VR,  A,  VN,  K,  J,  COMPV) 

This  subprogram  is  designed  to  compute  the  vector 
components  along  major  and  minor  axes  or  to  deter- 
mine north  and  east  components. 

DX,  VR — Angles,  magnitudes  of  popular  vector 

A — Azimuth  (A  =  0.  North  and  east  com- 
ponents) 

VN — New  component  series  which  is  formed 
K — Indicates  type  of  component  formed 
K  =  1.  Major  axis;  north  component 
K  =  2.  Minor  axis;  east  component 
J — Number  of  data  values  in  series 

COMPV — Magnetic  variation  in  directions  if  the 
directions  are  not  true  readings 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES 

ENVIRONMENTAL  SCIENCE  SERVICES  ADMINISTRATION  -  C.+  G.S. 

FOURIER  -  HARMONIC  ANALYSIS     PROJECT  13140?    29  DAYS 

PROGRAM  CUR  AN  ( I NPUT • OUTPUT • TAPF9 ) 

DIMENSION   VIN(9000)  ' DIN<9000 )  'CHSIpt ] 0)  'CRP( 10 ) ' COVU(25 ) 'CF( 25  >  ' 
1TXXIXU  5'20)'  7XXIX2  (  5  '  20  )  '  CK.APAP  (  20  )  '  CKAPI20)'  THS  I  (  20  )  '  CR  (  20  )  ' 
2PELL ( 20  )  ' BI NGO( 20 >  ' G( 5 )  ' PAK ( 5 >  ' D( 24)  ' W( 24 )  '  ZFTA(S)'  CRC ( 5 > 

DIMENSION  RES (4) ' ZEP( 4 ) 'XTEMP(9000) • VN(9000) ' ITEMPI9000) 

DIMENSION    CV(300)'     CD(300)'     ICop(300>'     IDBNSU8) 

DIMENSION  CXX (30)  ' DEG(40 )  '  TA&1(40)'  TAB2I40)'  TA63I40)'  TAB4(40) 

DIMENSION  TAB5(40)'  TAB6 ( 40 ) ' AO (4 ) '  BQ(3>'  OFX(^) 

DIMENSION  TXXIXA( 20)  ' TXXIXBI  20  >  ' TXXIXCt  20)  'TXXIXDI  ?0  > 'TXXIXEf  20) 

DIMENSION  TXXiXF(2  0)'TXXIXG(2  0),TXXIXH(2C),TXXIXl(7  0),TxXIXJ(20) 

COMMON / FOR T / AQ ' RQ/COSTX/ CXX' OEX/CONTT/NN 

COMMON  /  LEAPYR  /  DAY' MONTH 

SPEFDS  OF  CONSTITUENTS  M ( ? )  ' N ( 2  )  ' S ( 2  )  ' 0 ( 1  )  ' M  1  )   FOR  1  SAMPLE/HP. 

DATA  SXM2'SXS2' SXN2'SXOl' SXK1/0. 16 10228' 0.1666667' 0.1579985 "0.0774    15 

1613*0.08 3 5615/  16 


TABLE  XXIX  FROM  SCHUREMAN'  S.P.  °8    (194]) 

ACCELERATION  IN  EPOCH  OF  K(l)  DUE  TO  PC])'  20  DAYS 

DATA(TABl(JZ),J2sl'18)/11.4,16.4'18.3,17.6,I5.2,11.7,7.4,2.7'-2.2'      1 
1-6.9'  11.3,-14.9'-17.5,-18.3'-16.7'-12.1'-4.7'3.9/  2 

RATIO  OF  INCRFASE  IN  AMPLITUDE  OF  K(l>  DUE  TO  P(1)'29DAYS 

DATA(TAB2 (JZ)  ,JZsl,18)/-0.26,-0.17,-0.06,0.04,0.14♦0.23•0.28'0.3  1,      3 
10.32'  .29  •0.2?'0.15I 0.05 '-0.06 '-0.16 ' -0. 25 ' -0. 30 ' -0.31/  4 

ACCELERATION  IN  EPOCH  OF  S(2)  DUE  TO  K  (  2  )  '  29DAYS 

DATA!  TAB3  (  JZ)  '  JZsl,18)/5.9'Q.6,I2.6'l4.6,15.0']-;>.5*9.6,?.7,-3.01-9     5 
l.l,-13.2,-15.0l-14.7'-l2.9,-10.0'-6.4,-2.3'1.9/  6 

PATIO  OF  INCREASE  IN  AMPLITUDE  OF  S(2>  DUE  TO  K(2)'29DAYS 

DATA( TAG4( JZ ),JZsl'18)/0.24,0.19,0.]2,0.04,-0.05l-0.14'-0.21' 2*-0 .  7 
125' -0.21' -0.15 '-0.06 '0.03* 0.11' 0.18 '0.23 '2*0.26/  8 

ACCELFPATION  IN  EPOCH  OF  S(2)  DUE  TO  T<2)*  29  DAYS 

DATA( TAB5 (JZ) ,jZsl,37)/-0.8,-l.?,-1.8,-2.2'-2.6,-2.9,-3.2,3--3.3l-  9 
]3.]!,-2.8,-2«5'-?.0'-'|.5,-0.9,-0.3l0.3'0.9'1.5'2«0,2«4'2.8,3.1,?*3.     10 

23'^.2,3.0,2.7'2.^"l.'3,1.4,ri.8'n.3'_n.2,-C.8/  11 

RESULTANT  AMPLITUDE  IN  S(2)  DUE  To  T(2)'  29  DAYS 

DATA! TAB6( JZ)  ljZsl,37)/i.06,2*1.05'1.04,2*1.03,1.02,i.01,1.00,0.99  12 
l'0.98,0.97,0.96'2*0.95,4*0.94l2*0.95'0.9  6,0.97,0.98,0.99,l.O0'1.0l  13 
2 '1.02' 2*1. 03 '1.04' 2*1. 05' 4*1. 06/  14 


AMPLITUDE  EFFECT  OF  CONSTITUENTS  ON  Ml?)'  29  DAYS 

DATA(TXXIXA(  J  )  ' J s  1  '  20 ) / 0 . 0 ' 0  .  050 • 0 . 0 1 8 ' 2*0 . 0 ' 0 . 0 56 ' 0 .0 50 ' 0 . 049 ' 0  17 

1.021'  .060 '0. 096*0.0] 8'0. 096' 7*0.0/  18 

AMPLITUDE  EFFECT  OF  CONSTITUENTS  ON  N(2>'  29  DAYS 

DATA(TXXIXB(  J  )  ' Js 1 ' 20 > / 0 .0 50 ' 0 .00  '  0 . 00 5 ' 2*0 . 0 ' 0 . 052 ' 0 . 049 ' 0 . 050 '  19 

10. 031 '0.021 '0.018' 0.096' 0.968' 7*0.0/  2  0 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

AMPLITUDE  EFFFCT  OF  CONSTITUENTS  ON  S(2>'  29  DAYS 

DATA(  TXXIXC.I  J  )  '  Js  1  '  ?0  )  /0  .0 1  8  '  0  .005  '  3*0  .0  '  0  .  959  '  0  .096  '  0  .017  «  2*0  >  9     21 

190'0.  50' 0.'018 '0.042  '  7*0.0/  22 

AMPLITUDE  EFFECT  OF  CONSTITUENTS  ON  0(1)'  29  DAYS 

DATA(TXXIXD(  J  )  ' Js 1 ' 20 > /4*0  .0 ' 0 .0 56 ' 8*0  .0 ' 0 . 052 ' 0  .065 ' 0  .0 52 ' 0  .0 1 8     23 
1 ' 0.05  ' 0.049 ' 0.096/  24 

AMPLITUDE  EFFECT  OF  CONSTITUENTS  ON  Kf 1 ) '  29  DAYS 

DATA(TXXIXE(  J  )  ' J s ]  ' 20 ) / 3 *0 . 0 ' 0 . 0 56 ' 9*0  .  0 ' 2*0 . 050  *  0  .0 56 ' 0 . 9 59 ' 0 . 0  25 

152'0.49'O.Ol]/  26 
PHASE  EFFECT  nF  CONSTITUENTS  ON  Mfp)'  29  DAYS 

DATA(TXXIXF(  J  )  '  Js  1  ' 20 ) /0  .0 ' 3 5  1  . ' 1 74  .  '  2*0  .  0  '  22  .  ' 9  .  '  34  1  .  '  8  .  '  1  59 . ' 1  27 

164. '  186.  •  196.  ' 7*0.0/  28 
PHASE  EFFtCT  Of  CONSTITUENTS  ON  N(2)'  29  DAYS 

DATA(TXXIXG(  J  )' Js 1 ' 20 ) /9 .' 0  .  0 ' 3 .' 2*0  .  0  '  32 .'  19  .'  35  1  •'  1  7  .'  1  69 .'  1 74  29 

1 . ' ] 96. ' 25. ' 7*0.0/  30 
PHASE  EFFECT  OF  CONSTITUENTS  ON  S<2>'29  DAYS 

DATA(TXXIXH)  J  )  ' Js 1 ' 20 ) / 1 86 . ' 357 . ' 3*0 .0 ' 29 . ' 196 . ' 348 . ' 14. ' 346  .  ' 3 5  31 

11.' 193. '202. ' 7*0.0/  ^2 
PHASE  EFFECT  OF  CONSTITUENTS  ON  0(1)'  29  DAYS 

PATA(TXXIXI(  J  )  ' Jsl ' 20) /4*0.0 ' 22. ' 8*0.0  '  32. ' 13. '44.  ' 174. ' 351  .  • 341  33 

] . ' ] 96./  34 
PHASE  FFFFCT  OF  CONSTITUENTS  ON  K(l)'  29  DAYS 

DATA(TXXIXJ(  J  )  ' Jsl ' 20 ) /3*0.0 ' 338. ' 9*0.0 '9.  '  351.  ' 22.  ' 331. ' 328.  ' 31  35 

19. '3  54./  3  6 

LOADING  OF   TABLES  INTO  PPOPEP  ARRAYS 

DO  80  8  MAXJ  s  1 ' 20 

TXX  1X1  (  1  'MAXJ  )  s  TXXIXA(MAXJ) 

TXXlXl  (  ?'MAXJ  )  s  TXXIXP(MAXJ) 

TXXIX1  (  ?'MAXJ  )  s  TXXIXT(MAXJ) 

TXXIX1 (4'MAXJ)  s  TXXIXMMAXJ) 

TXX  IX]  ( 5'MAXJ  )  s  TXXIXF(MAXJ) 

TXXIX21  1  'MAXJ  )  s  TXXIXF(MAXJ) 

TXXIX2 ( 2'MAXJ )  s  T XX  I  XG ( M AX J ) 

TXXIX2I 3'MAXJ )  s  TXXIXH(MAXJ) 

TXXIX2 (4'MAXJ )  s  TXXIXI(MAXJ) 

8008  TXXI  X2( 5'MAXJ )  s  TXXIXJ(MAXJ) 

CONTROL  CARDS  ENTERED 

READ  221 'NJ 

DO  1111  JN  s  1'NJ 

READ  808' ( IDENSI JK) ' JK  s  1'18> 

READ  76'  AZI ' LSORT1 Jj'N'MAP' K ' CVAR ' INC 'NSPH' TFAC'EAZI 

RFAD  2038' XYEP  'MONTH' DAY ' STT ' TM ' GONL 

RFAD  2039'  AQ '  BO 

IF ( NSPH.EQ.O.OP.TFAC.EO.O.O)  GO  TO  8356 

DATA  INPUT  LOOP  BEGINS 

GO  TO  (  101  '  102  '  103' 660 • 661 )' MAP 
661  READ  AQ,' (DIN(KX) 'VIN(KX)  'KX  s  1'N) 

GO  TO  104 
101  READ  AQ1 ( VINI I ) ' I  s  1 'N) 

GO  TO  123 
103  PFAD  AQ'  ( VIN ( I ) '  I  s  1'N) 

IF ( AZI .EO.O.O.AND.K.FQ.O )  GO  TO  123 

READ  BQ'  (DIN<  I  )  •  I  s  1'NI 

GO  TO  104 
123  IF(LSORT  -  1)  120'12]'l21 
120  CALL  TEMPI VIN'XTFMP'N) 

GO  To  108 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

1?1  CALL  SORT  I VIN'XTFMP'N' JJ ) 

N  s  NN 

GO  TO  108 
10?  READ(9'AQ)  (DIN(I  )  '  VIN(I)'I  s  l'N) 

IF(K.EQ.5)  GO  TO  120 

GO  TO  104 
660  READ  222'MX ' LPD' LX 

IF(LX.FO.O)  GO  TO  666 

RFAD  555 '  (  ICOR (LY)  'Cp(LY)  'CV( LY  )  ' LY  s  1  ' LX ) 
666  CALL  EDAT(MX' VIN'DIN'N'LX' ICOR'CD-' CV'LRD'XTFMPi  INC) 
104  NFCOM  s  1 

IF(K.FQ,3)  NECOM  s  3 

IF(K.FQ.4)  NECOM  s  4 

IF(MECOM.EG.3.0R.NECOM.E0.4)  K  s  1 

IF ( NECOM. EQ. 4 )  GO  TO  8155 

CALL  DETAZI DIN' AZI ' N • CvAR •  VIN ' XTEMP'  INC) 

FLX  s  AZI 

IFIFAZI.  EQ.0.0)  GO  TO  8355 

CALL  DETAZI DIN' EAZI 'N'CVAR'VIN'XTEMP'O) 

FPX  s  FAZI 

DIPF  s  FLX  -  FBX 

AT7  s   .0 

PIPF  s  APS (DIPF) 

I  F  I  P  I  P  E  .  G  T  . 180.0)  ATZ  s  PIPF  -  180.0 

IF( PIPE. LT.l 80.0)  ATZ  s  180.0  -  PIPF 

APZ  s  ATZ/2.0 

ZAP  s  AINT(ARZ) 

IF (DIRE. LT. 0.0  )  AZI  s  ELX  +  ZAP 

IF(DIRE.GT.O.O)  AZI  s  FLX  -  ZAP 

PRINT  7359'  AZ I ' DI R E • AR Z ' Z AP 

8355  CALL  AZIM  (  DIN  •  VIM  '  AZ  I  '  VM  '  K  'N  '  CVAR  '  NFCOM  ) 
I  F ( LSORT  -  1 ) 107 ' 106 ' 106 

106  rALL  SORT ( VN ' XTfmp' N ' JJ) 
N  s  N  M 

GO  TO  108 

107  CALL  TEMP (VN ' XTFMP'  M) 

FOURIER  ANALYSIS  PERFORMED  FOR  FIVF  MAJOR  TIDAL  FREQUENCIES 

108  XM?  s  TFAC*SXM2/FLOATf NSPH) 
XS?  s  TFAC*SXS2/FLOAT ( MSPH) 
XN?  s  TFAC*SXN2/FLOAT (NSPH) 
XO  s  TFAC*SX01 /FLOAT ( NQPH ) 
XK  s  TFAOSXSO  /FLOAT  (N^PH) 
GO  To  8^57 

8356  PRINT  7358 
CALL  FXIT 

83^7  CALL  FOPAN ( XTEMP ' XM2 ' RFS ' ZEP ' N  '  4  ) 
CHSIPI  1  )  s  ZEP(1  ) 
CHSIP(6)  s  ZEP12  ) 
CHSIPI 7  )  s  ZEP( 3) 
CHSIP(8)  s  ZFPI-4) 
CPP (  1  )  s  PES (  1  ) 
CRD ( m  s  RES ( 2  ) 
CPP(7)  s  R  E  S ( 3 ) 
CPP(8)  s  PES (4) 

110  CALL  FORAN(XTFMP ' XN2 '  PES'ZFP'  N'  1) 
CHSIPf  2)s  ZFP (  1  ) 

CPP ( 2  )  s  RFS (  1  ) 

111  CALL  FORAN(XTEMP'XS2'  RES'ZFP'  N'  3) 
CHSIPI 3 )  s  ZFP( 1 ) 

CHSIP(Q)  s  ZEP( 2) 

CHSIPI 10)  sZEP( 3) 

CPP( 1)  s  RES! 1 ) 

CRP( 9)  s  PES( 2  ) 

CPP (lis  RFS ( 3) 

CALL  FORAN(  XTEMP'XO"  PFS'  ZFP'  M'  1) 

CHSIP(4)  s  ZFP(1 ) 

CPP (4)  s  RES( 1 ) 

CALL   FORANtXTFMp ' XK'  PES'  ZFP'  N'  1) 

CHS  IP! 5 )  s  ZEPI 1  ) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES -Continued 


CRP(5)  s   RES(l) 
CALL  TWOPI (CHSIP' 10) 

IF(NEC0M.EQ.3.0P.NECOM.E0.4) 
GO  TO  8057 
R05R  IF(K.FQ.2)  GO  TO  8059 


GO  TO  8058 


DETERMINATION  OF  VO+  U  FOLLOWS  WITH  LOGF  AND  ARGUMENTS 


R057  CALL  DAYXXtMONTH' DAY'STT" TM' DAyB'DAyM' GRES'GPMS'CP 

DO  7135  JAN  s  1'19 

LXN  s  JAN  +  18 

TAP!  ,'LXN)  s  TAB!  (JAN) 

TAR2(LXN)  s  TAB?(JAN) 

TAB^(LXN)  s  TAB3(JAN) 
71^5  TA«4(LXN)  s  TAB4UAN) 

CALL  ASTPOIXYBR'DAYB' DAYM • GPBS ' GRM5 • 29) 

PM  s  CXX( 1  ) 

PL  s  CXX(2) 

<^L  s  CXX(3) 

PS  s  CXX(4) 

PLM  s  CXX(5) 

SKYN  s  CXX(6) 

VI  s  CXX19) 

V  s  CXX( 1 0  ) 

XI  s  CXX( ]  1  ) 

VP  s  CXX( 1 2  ) 

VPP  s  CXX ( 1 3 ) 

P  s  0XXU4) 

AUL  s  CXX115) 

AUM  s  CXX (16) 

PAPA  s  0.0 

DO  3073  NAP  s  1 '37 

DFG(NAP)  s  PAPA 
^073  PAPA  s  PAPA  +  10.0 


DETERMINATION   OF  EQUILIBRIUM   ARGUMENTS  (VO  +  U) 


7777  TML  s  CP  -  GOML 
COM  s  SL  +  TML 
COVU( 1 )  s2.0* (CON  - 


COVUt?) 

covU(3 ) 

COVU(4) 
COVU(5' 
C0VU(6) 
COVLU7) 
COvU(8 ) 
COVU(9 ) 
COVUt 1 0 
COVU( n 
COVU( 1 2 
COVU( 13 
COVlM  14 
COVU( 15 
COVU(16 
COVU( 17 
COVUt 18 
COVUt 19 
^OVU(2C 
COVUt 21 
COVUt 22 
COVUt 23 
COVUt  24 
COVU( 2  5 


PM  +XI  -V ) 
PL) 


S  COVUt 1 )  -( pm  - 

s  2.0#(TML) 

s  CON  -  V  -  2.0* (PM  -  XI)  -  90.0 

s  CON  -  VP  +  90.0 

s  2.0 -"CON  -  VPP 

S  2.0*  CON  -  PM  -PL  +  AUL 

s  COVU(?)  -(pm  -  PL) 

S  SL  -  Pc,  +  180.0  +  2.0*  (TML) 

s  COVUC*)  -  (SL  -  PS) 

s  COVLM  i)  -(2.0*(SL  -  PM))  -  ( PM  -  PL)  +  180.0 

s  COVLM 1 )  +  2.0* ( SL  -  PM) 

s  COVLMl2>  +  (PM  -  PL) 

s  CON  +  PM  -  PL  -  V  +  90.0 

S  CON  -  PM  +  AUM 

s  CON  -  V  +  2.0* (PM  -  XI)  +  90. C 

s  TML  +  270.0  -  SL 

s  C0VU(4)  -  (PM  -  PL> 

s  COVIM18)  -  (PM  -  PL) 

s  COVU( 18>  +  ?.0*( SL  -  PL ) 

s  2.0* (COVLM 1  )  ) 

s  3.0*  (COVLM  1  )  ) 

s  4.C*(C0VU( 1 ) ) 

s  2.  0*(<~OVU(  3  )  ) 

s  3.0*(COVU( 3) ) 


CALL  TWOPI  (COV'J1  2^  > 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

DETERMINATION   OF  NODE  FACTOR  RECIPROCALS  (F) 

V  s  V*. 0174533 
P  s  P*.  01 745 3 3 

VI  sVI -'.0174533 

CPA  s  1.0/SQRT(1.0  -  ]?.0*(SIN(.5*VI )**2/COS(.5*VI )**2)*COS(?.0*P) 
1+  if..     *(  SIN  (  .S"VI  )  »*4/00S  (  ,SJ"-VI  )**4  )  ) 

COA  s  1.0/SQRT(?.3lC  +  ] . 43 5* COS ( 2 . 0*P ) ) 

CF(1)   S  1  .0/  (  r^c  (  .5*VI  )"*4/.91  ^44  ) 

CF(?)  s  CF(1) 

C  F ( 3 )  s  1.0 

C|r(4)  s  1  .0/  ('SIM  (VI  )*COS(  .5*VI  )**2/.  37988  ) 

CF(5)  s  1  ,0/SORT  (  0.8°65"-c  INI  2.0*VI  )  **2  +  0  .  6001*SI  N  (  2  .0*VI  )  *COS  (  V  ) 
1+0.1   6) 

CF(6)  s  1.0/SQRTI 19. 0444*51 N < VI  )**4  +  2  .  77C2*SI N ( VI ) **2*CCS ( 2 .0*V ) 
1+  0.0981)  -  0.000CJ99 

CF (7  )  s  CF( 1 ) *CRA 

CF(8)  s  CF(1) 

CF(9)  si .0 

CF( 10) s  1.0 

CFIllls  CF(1) 

CF(12)  s  rF  n  ) 

CF( 13)  s  CF(1 ) 

CFU4)  s  1.0/ (SIN(2«0*Vl ) /0. 72137) 

CF ( ] 5 )  s  CF (4 )*COA 

CF(16)  s  1.0/(SlN(VI)*SIN( .5*VI )**2/0. 016358) 

CF(  17)  s  1.0 

CF(  18)  s  CF(4) 

CF{  19)  s  CF ( 18  ) 

CF  (  20)  s  CF( 19) 

CFI21)  s  CF(1)**2 

CF(22)  s  CF( 1 )**3 

CF(23)  s  CF ( 1 )**4 

CF(24)  s  1.0 

CF(25)  s  1.0 

CKKF  s  1.0/CFI6) 

AHNAT  s  CF(4)*CRP(4) 

C5F  s  CF(5  ) 

ENTER   SECONDARY  EFFFCTS  AND   DETERMINE  KAPPA 

CALL  SASTRfDEG'CjF'CKKF' SMAC 'PROD' ACCP' RESAM « 37 ' TAB1 « TAB2 ' TAB3 ' TAB 
14,TAR5'TAB6,SL'VPP'Pc,VP'OF'?5) 
8059  DO  21  I  s  1  '  c 

2]  CKAPAP(I)  s  CHSIP(I)  +  COVU(I) 
CALL  TWOP I  ( CKAPAP'  5  ) 
CKAPt] '  s  CKAPAP ( 1 ) 
CKAP( 2  )  s  CKAPAP ( 2  ) 
CKAP(3>  s  CKAPAP!  3)  +  SMAC 
CKAP(4)  s  CKAPAD(4) 
CKAP(5>  s  CKAPAP(5'  +  ACCP 
CKAP (6 )  S  CKAP ( 3 ) 
DMN  s  CKAP( 1 ) -  CAP ( ? ) 
CALL  ONEP I (  DMN) 
CKAP  (  7  )  s  CKAP  (  1  )+  D«N 
CKAP(8>  s  CKAPt?)-  DMN 
CKAP (9>  s  CKAP (3) 
CKAP( 10  )  s  CKAPf  3 ) 
SMD  s  CKAP( 3) -  CKAP( 1 ) 
CALL  ONEPI (SMD) 

CKAP(ll)  s  CKAP(l)  +  0.464*SMD 
CKAP(12)  s  CKAP( 1 )  -  SMD 
rKAP(13)  s  CKAP(l)  -  0.866*DWN 
OKn  s  CKAPi5)  -  CKAP (4) 
CALL  ONFPI (OKD) 
CKAPI14)  s  CKAP(5)  +  0.5*OKD 
CKAPM5)  s  CKAP(5'  -  0.5*OKD 
CKAP(  16  )  s  CKAP( 5)  +  OKD 
CKAP(17)s  CKAPI5) 
CKAPU8)  s  CKAP(5)  -  1.5*OKD 
CKAP(19)  s  CKAPI5)  -  2.0*OKD 
CKAPI2C)  s  CKAPI5)  -  1.43*OKD 
CALL  TWOP I (CKAP' 20 ) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

DETERMINATION  OF  ZETA(A)  AND  R(A)  FOLLOWS 


DO  2  2  I  s  1 ' 20 
22  CHSI(I)  s  CKAP(I)  -  C0VU'(  I  ) 
CALL  TWOPI ( CHSI •  20) 
CR  (  1  )  s  CRP( 1 ) 
CR (2  )  s  CRP( 2  ) 
CR ( 3 )  s  CRP( 3 ) /PPOD 
CR (4)  s  CPP(4) 
CR(5)  s  CRP(5)/RFSAM 
CR(6)  s  0.272*CP ( 3  ) /rF(6  ) 
CP(7)  s  0.145*CP ( 2 ) /ORA 
CR(8)  s  0.133*CR  (  2  > 
CP(o)  s  0.008*CP<3) 
CR ( 1 0) s  0.059*C? ( 3 ) 
CP( ] ] )  s  0.007*CP( 1 ) 
CR(12)  s  0.024*CR(1) 
CR(  1?)  s  0.194*CP(2 ) 
CR(14)  s  0.079*AHNAT/CF( 14) 
CR(15)  s  0.071-*CR(4)/CQA 
CP(16)  s  0.043*AHNAT/CF < 16) 
CR (  17)  s  0.331*CP ( 5)*CF(  5 ) 
CR (] 8  )  s  0.194*CP ( 4 ) 
CR (  19)  s  0.026*CP ( 4 ) 
CP(20)  s  0.038*CP(4) 


ELIMINATION  OF  COMPONENT  EFFECTS  AND  DETERMINE  KAPPA  AND  H(A: 


DO  35  I  s  1' 5 
SUM  s   .0 
SOO  s   .0 
DO  20  J  s  1'20 
311  BELL(J)  s  CR ( J)*TXXIX1  (  I  '  J) 
IF( BELL (J)-. 0005)  23' 24 '24 
23  BELL(J>  s  0.0 

iu    BINGO(J)s  TXXIX?tI'J>  -  CHSI(J)  +  CHSIP(I) 
IF (BINGO (J) )  2^' 2ft ' 26 

25  BINGO(J)  s  <  RINGO(J)  +  360 .) *0 . 01 745^3 
GO  TO  29 

26  IF(BINGO(J)  -  360.0)  28'28'27 

27  BINGn(J)  s  (BINGO(J)  +  360.0  1*0.0174533 
GO  To  29 

28  BINGO(J)  s  BI NGO( J  1*0.0174533 

29  SUM  s  SUM  +  BELL ( J )*SIN( BINGOI J) ) 
?0  SOO  s  SOO  +  BELL ( J>*COS( BINGO( J) ) 

QOS  s  CPP( I )  -  SOO 

CALL  TI TAN (SUM 'OOS1 DCCH' 2 ) 

CDCH  s  DCCH*0 .017453^ 

^pc  (  i )  s  o.o 

IF ( DCCH. EQ. 90.0. OR. DCCH. EO. 270.0)  Go  TO  34 
CRC(I)  sf^RPlI)  -  SOO) /COS(CDCH) 

34  G(  I  )  s  CRC(  I  ) *CF(  I  ) 

ZETA(I)  s  CHSIP(I)  +   CDCH*57. 29578 

35  PAM  I  )  s  CHSIP(I)  +  CDCH*57. 29578  +  COVU(I) 
CALL  TWOPI ( PAK '  5 ) 


INFEPFNCE   OF       H(A)      FPOM  MAJOR  CONSTITUENTS 


D(  1  )  s  0.079*G (4 ) 

D  (  7  )     s  0(5) 

D(3)  s  0.272*0(3) 

DIM  s  0.-O2  8*G  (  1  ) 

D  (  5  )  s  0  .  0  7  ]  *  G  (  4  ) 

D  (  6  )  s  G  (  ]  ) 

D(7)  s  CPP (6 ) *CF ( 21 ) 

D(8)  s  CPP ( 7 ) *CF ( 22 ) 

D(Q  )  s  CRP (8  )  *CF( 23 ) 

niiO)  s  G  (  2  ) 

D<  1  1  )  s  0.133*G(?J 

n  ( i  2  )  s  G  (  4  ) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 


D(l^)  s  0.043*G(4) 

Dt  14)  s  0.331*G< 5 ) 

D(15)  s  0.194*G(4) 

D(16)  s  0.026*G(4) 

D(17)  s  0.008*G( 3 ) 

D  (  ]  8  )  s  G  (  3  ) 

D<  1  °  )  s  CRPI 9  )  *CF< 2  4 ) 

D(?0)  s  CRP ( 10)*CF(2S) 

DC?]  )  s  n.059*G( ? ) 

D(?2 )  s  0.007*01 1 ) 

D(?3)  s  0.194*G(?) 

D(?4)  s  0.038*G(4) 


INFERENCE   OF 


KAPPA(A)    FROM.  MAJOR 


)NSTITUENTS 


677<v 


6780 
6781 


DAFX 

GAFX 

PAOT 

FPAX 

PATO 

PATTO 

CALL 

CALL 

IF(AE 

I  F  (  GA 

FPAX 

GO  TO 

FPAX 

CALL 

CALL 

W(1  1 

W  (  7  ) 

wo 

W(4) 
W(  "5  ) 
W(6) 

W  (  7  ) 
W(8) 
W  (  o ) 
W(  1  0 
W  ( 1  ] 
W(]  2 

wen 

W<  1  4 
W(  15 
W<16 
W(  17 

wn  8 

W(19 

Wt  2  0 
W  (  7  1 
W(?2 
WC?3 
W(  24 
CALL 


PAK ( 5  ) 

PAK ( 5  ) 

PAK (  1  ) 

PAK ( 5  ) 

PAK  !  2  ) 

s  PAK ( 3 

ONFP  I  ( PA 

ONEP  I  ( PA 

S(GAFX) . 

FX  16779' 

s  FPAX  + 

678] 

s  FPAX  - 

ONEPI ( PA 

ONEP I ( PA 

s  PAK (5 

s  PAK ( 5  ) 

PAK ( 3  ) 

2.0*PA 

0.5*FP 

PAK(  1  ) 

CHSIPI 

CHSIPt 

CHSIPf 

PAK  (  2 

2.0*P 

PAK(4 

2.0*P 

PAK  (  5 

PAK(  5 

PAK  (  5 

PAK(  3 

PAK  (  3 

CHSIP 

CHS  IP 

PAK  (  3 

PAK  (  3 

PAK  (  1 

PAK  (  5 

TWOP I ( W 


)  - 
FX) 
OT) 
LT. 

678 
36 


PAK (4) 
PAK (4  ) 
PAK  ( 2  ) 
PAK U  ) 
PAK  (  1  ) 
PAK  (  1  ) 


180.0)  GO 
1  '  6  7  8  n 
0.0 


TO  6781 


360.0 
TO) 
TTO) 
)  +  0.5*PAFX 


K ( 1  )  -  PAK ( 2  ) 
AX 


+  COV'.'(21  ) 

+  COVU( 22  ) 
+  C0VUI23) 

(?)     -    PAK  (  1  ) 

(5)  -  PAK (4) 

-1.5*PACX 
-  2.C*PAFX 


COVU( 24  ) 
COVLM  25) 


-  0.536*PATTO 

-  0.866*PAOT 

-  1.43*PAFX 
) 


STATION  IDENTIFICATION  PPINTFD 


PRINT  13 

PRINT  810'  (  I  PENS! K J)  ' KJ  s  1 ' IS ' 


ZFTA  PRIME  AND  P  PRIME  OUTPUT  FOR  FIVF  MAJOR  TIDAL  CONSTITUENT! 
AND  SECOND' THIPD1  AND  FOUPTH  HARMONICS 


PRINT 

811 

PR  I  NT 

73' 

(CHSIPf I ) 

'I  s  1  '  5) 

PRINT 

812 

(CRPI I ) ' I 

s  1  '  5> 

PRINT 

813 

PRINT 

73' 

(CHSIP (  I  ) 

'I  s  6' 10) 

PRINT 

812 

(CRPI  I  )  ' I 

s  6' 10) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 


PRINT  5 

PRINT  6'  CKKF 'CPA'CQA  '  AHNAT 


EQUILIBRIUM  ApGU^ENTS  AND  NODE  FACTOPS  DnINTED 


PPIMT 
PRINT 
PRINT 
PRINT 
PRINT 
PPI  NT 
PR  I  NT 
PPI  NT 
PPI  NT 
PRINT 
PRINT 
PRI  NT 
PRINT 
PRINT 
PRINT 
PRINT 


7 

811 

814 

815 

816 

814 

815 

817 

814 

815 

818 

814 

815 

813 

814 

815 


!COVU( I 
(CF(  I  )  ' 

(COVU( I 
tCF(I)  • 

(covut  i 

(CF(  I  )  ' 

<COVU(  I 
(CF(  I  )  ' 

(COVU( I 
(CF (  I  )  ' 


5  ) 

6'  II 

10) 


I  sll'15) 
s] 1  '  15) 

I  sl6' 20) 
S161 20  ) 

I  521*25) 
s21 '25) 


KAPPA"  ZETA'  AND  P  PRINTED  FOR  TIDAL  CONSTITUENTS  AND  HAPMONICS 
PRIOP  TO  ELIMINATION  PROCESS 


PRINT 

PRINT 
PRINT 
PPI  NT 
PPI  NT 
PRI  NT 
PPI  NT 
PRINT 
PRINT 
PRINT 
PPIMT 
PRINT 
PPI  NT 
PRINT 
PRI  NT 
PRI  NT 
PRINT 


9 

811 

819 

820 

821 

816 

819 

8  20 

821 

817 

819 

820 

821 

818 

819 

820 

821 


( CKAP( I ) 
(CHS  I  (  I  ) 
(CR  (  I  )  •  I 

(CKAP( I ) 

( CHC  1(1) 
(CR  (  I  )  '  I 

(CKAP( I  ) 

(CHS1  (  I  ) 
(CR (  I  )  •  I 

(CKAP( I ) 
(CHSI ( I ) 
(CR  (  I  )  '  I 


I  s  1 '  5) 
I  s  1 '  5) 
s  1  •  5) 

I  s  6' 10) 
I  s  6'  10) 
s  6'10) 

I  sll'15) 
I  sll'15) 
sll'15) 

I  sl6'20) 
I  s!6'20) 
sl6'20) 


R'  ZETA'  KAPPA'  AMD  H  PRINTED  FOR  FIVE  MAJOR  TIDAL  CONST  I TUENT S 
AFTED  ELIMINATION  PRocfSS 


PPINT 

11 

PRINT 

822' 

PRINT 

823' 

PRINT 

824' 

PRINT 

825' 

PPIMT 

826' 

CRC(  1  )  '  ZETAI  1  )  '  PA!:"  (  ]  )  '  G  (  1  ) 
CPC(2)'ZETA(2)'PAK(2»'G(2) 
CRC13)  'ZETAI 3)  ' PAK ( 3>  'G(3) 
CRC(4)'ZETA(4)'PAK(4)'GU) 
CRC(5)'ZETA(5)'PAK(5),G(5) 


COMPUTED  VALUES  nF  H(A)  AND  KAPP.A(A)  nPINTED  FOP  THE  CONSTITUENTS 

M2    S2  N2  01    Kl 

M4    M6  M8  S4     S6 

INFERRED  VALUES  OF  H(A)  AMD  KAPPAIA)  PRINTED  FOP  THE  CONSTITUENTS 

K2    L2  2N  R2    T2    LAMPDA    MU2    NU2 

Jl    Ml  001  PI    01    ?01    PH01 


PRINT 

15 

PPINT 

16 

PPIMT 

71  ' 

(  W  I  I  )  ' 

I  s  1 

5) 

PPIMT 

72' 

( D(  I  )  ' 

I  s  1 

5) 

PPINT 

17 

PPINT 

7] 

(W(  I  )  ' 

I  s  6 

10) 

PRINT 

72 

( D(  I  )  ' 

I  s  6 

10) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 


5899 
1112 
7111 


PP  I  NT 

18 

PRI  NT 

71 '  (W(  I  )  ' 

I  s 

n 

1  5  ) 

PPI  NT 

72'  ( D(  I  )  ' 

I  s 

11 

1  5) 

PRINT 

19 

PPI  NT 

7  3  '  (  W  (  I  )  ' 

I  s 

16 

20) 

PRINT 

7  2-'  (  D(  I  )  ' 

I  s 

16 

?0) 

PRINT 

220 

PRINT 

7  1  '  (  W  (  I  )  ' 

I  s 

21 

24) 

PPI  NT 

7360 

PPI  NT 

72' <D( I ) • 

I  s 

21 

■24) 

CALL  XMEANIXTFMP'SUTZ'N) 

IF  (  NFCnM.F0.3  .AMD.K.FO.  1  )  GO  TO  5°°R 
IF(NFCCM.EQ.4.AND.K.FQ.1>  GO  TO  5998 


GO  TO  5899 
5998  K  s  2 

GO  TO  (5899'5 

<S9Q9  AZI  s  90.0 

GO  Tn  8355 

GO  TO  (1111 

REWIND  9 


S99'  8355'  0355  )  '  NEC 


1112'1111'  1112'  1111  )  'MAP 


CONTINUE 
STop 


5 
6 
7 
9 
11 

13 

14 
15 
16 
17 

1  R 

10 

220 

71 

72 

73 

74 

76 

77 

78 

79 

221 

222 

5^5 

808 

81  G 

81  1 

812 

81  3 

87  4 

815 

81  6 

817 

878 

879 

820 

821 

822 

823 

824 

825 

826 

?0?7 

2038 

2039 

4O4O 

4041 

4440 

7358 
7359 
7360 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
14H7FTA 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FODMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FOPMAT 
FORMAT 
END 


//  28X'7H1/F( 

1H  '25X'  F9.4 
21X'49HEQ'JILI 
/  42X' 22HANAL 
1H0'3X'33H  EL 
5X'  -5HKAPPA' 
1H1  '48X' 27HHA 
1H1'4CX'27HHA 
///  14X'19HHA 
16X'4HJ(1 ) '6X 
16X*4HM(2) ' 6X 
16X'4H ( ?M)  ' 6X 
]6X'4H(2Q! '  6X 
16X'4HT (2 ) '4X 
5X'  5HK.APPA' 
6X  '  4HH( A )  • 
9X' 11HZETA1P 
24F3.1  ) 
F5.0'5I5'F5.0 
6CX'  F3.C  24X' 
24F3.2  ) 
24F3  .0  ) 
15) 
315) 

15'FlO.O' FlO. 
9A8/9A8  ) 
/24X'9A8//24X 
27X'4HM( 2  )  ' 8X 
1H0' 10X  '  RHR  (P 
27X'4HM(4) ' 8X 
1  2  X  '  7  H  (  V  +  U 
1H0 ' 13X '  4HF  ( 
27X '4Hk ( ?) ' 8X 
6HLAMBDA' 
5HOOI7  ) '8 
15X  '  5HKAPPA' 
1H0M5X'  4HZE 
15X1  4HR( 
4HM( 2 ) 
4HN( 2 ) 
4HS( 2 ) 
4H0( 7 ) 
4HK ( 1 ) 
5F7.3  ''F5.0  ) 

F5.0'  I  5* F5.0 
4A8'  3A8) 
1  HO '  lCX  '  F6.0 
1H0'  4r10.5> 
1H0'  1  OF  12. 4 
1H0' 23HJOB  'TE 
5X'  F10.1  '4X' 
1H0) 


K2>  ' 

'  2F1 
BP  I  UM 
Y5I5 
IMINA 
5X'  4 
RMONI 
RMONI 
RMONI 
'4HK  ( 
•4HM( 
1  4110  ( 
•4HR( 
'  6HLA 
5^10 
5  F 1 0  . 
RIME  ) 


•  2HP 
5'  F 
V  + 
D  'IN 
ON  C 
(A) 
ANAL 
ANAL 
CONS 

•  6X' 
•6X  ' 

•  6X' 
'6X- 

da1  5 

/) 

//) 

5F12 


A'  8X '  2HOA'  AX •  6HHhO( 1 )  ) 

10.4//) 

U)  AND  ELIMINATION  (F)  ARGUMENTS  /) 

FERENCF  //) 

F  COMPONENT  rFFFCTS  // 1 5X ' 4HR ( A  )  • 5X 

) 

YSIS'   29-DAY.S  ) 

YSIS'  15-DAYS'  2X'13HN(2)  INFERRED) 

TANTS '  IX '  13H(H)  AND  KAPPA  //) 

4HK ( 2 )  '  6X ' 4HL ( 2 )  ' 6X  *4HM ( 1  )  ) 

4HM(6 ) • 6X'4HM(8 ) '6X*4HN(2 ) ) 

4H(0O) ' 6X'4HP(1 ) '  6X'4HQ< 1 ) ) 

4HS(2)  '6X'4HS(4) '6X*4HS(6)  ) 

X'  5HN'J  (  2  )  '  4X'6HRI-IO(  1  )  ) 


.3) 


'215' F10.8'F5.0) 
F4.2) 


2) 


26X1 
26X1 


1H0  ' 
5X1 
5X1 
5X1 
5X1 
5X' 


•  9A8///) 

'  4HN( 2 )  ' 8X'4HS(2)  ' 8X'4HO( 1 )  '8X'4HK (1  )  ) 

PI  ME)'  SF1?.3//) 

' 4HM( 6)  ' 8X  '  4HM( 8 )  ' 8X '4HS (4)  ' 8X'4HS (6 )  ) 

)  '  IX  •  5F12.3) 

A)  '2X'  5F12.5  // ) 

' 4HL(  2)  ■ 8X  '4H( 2M)  '8X'4HR(2)  '  8X'4HT ( ?  )  ) 

6X  '  5HMU ( 2  )  '  7X '  5HNU ( 2  )  '  8X ' 4hJ (  ]  )  ' 8X ' 4HM (  1  )  ) 


X '4HP( 1 ) ' 8X '4HO( 1 ) 
5F12.3) 


RX '4H( 2Q)  ' 7X '6HRHO (  1  )  ) 


TA" 

A)  ' 
FlO 

FlO. 5' 
FlO. 5' 
F70.5  ' 
F70.5' 


5F12.  1) 
5F12.4  //  ) 
5'  2F10.2' 
5'  2F10.2' 
5'  2F10.2' 
5'  2F10.2' 
2F10.2 ' 


FlO. 5) 
FlO. 5) 
FlO. 5) 
FlO. 5) 
FlO. 5) 


' F5.2' 2FA.2' 

'  6F10.2) 

///  IX'  10F12.4) 

RMINATED '24H  NSPH  OR  TFAC  s  0  (ZERO)) 

Fl0.1'4X'  2F10.1) 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SUBROUTINE  FOPANISER'  SP«  PES'  ZFP'  N «"  M  ) 

FOUPIER  ANALYSIS  -   DETERMINE  ZETA (PR  I  ME ) AND  O(PPI^E) 

DIMENSION   SER(N)'  RFS(M)'  ZFP(M) 
DO  55  JR  s  1'  M 
PES(JR>  s  0.0 
ZFP(JP)  s  0.0 
C  s  0.0 
S  s  0.0 

DO  10   L  s  1*  N 

C  s  C    COS(FLOAT(L-l)*<3.141593*(FLOAT<JR)  *SP  )  )  )  *SER  (  L  ) 
100  S  s  5    SIN(FLOAT(L-l )  « ( 3 . 1 4 1 593* ( FLOAT ( JP)  *SP ) ) ) *SER I L ) 
C  s  2.0/  FLOAT(N)*  C 
S  s  2.0/  FLOAT (M)*  S 
X  s  s/c 
IF  DIVIDE  CHFCK  14'15 

14  IF(S)16'55' 17 

16  ZEP( JR) s  270.0 
GO  TO  18 

17  ZEPI  JP )  s  90.0 
GO  TO  18 

15  ZFP(JR)  s  ATAN(X)*57. 29578 
IF(S)  32' 32' 34 

32  IF(C)  33'55'37 

33  ZED(JP)  s  180.0  +  ZFD(jR) 
GO  To  ]8 

■*4  IF(C)  35'17'36 

^5  ZFP(JP)  s  ZEP(JP)  +  180.0 

GO  Tn  18 
"*7  ZFP(JP)  s  ZEP(JR)  +  360.0 

GO  TO  18 
36  ZFP(JP)  s  ZEP(JP)*1.0 
18  PES(JR)  s  SORT!  C**2  +  S**2> 
55  CONTINUE 

PFTIIPN 

FND 


5URROUTINE  SASTP(DEG  C5 F ' CKK F ' ACS ' APS ' ACK ' Apr •  I OT ' T ABl  ' T AB2  '  TAB3  ' 
1TAR4,TAB5'TAB6'SL'VPP,PS,VP,<'F,JOP) 

DETERMINE  SECONDARY  CONSTITUENT  EFFECTS 

DIMENSION  DEG(I0T),TaBi(IOT)'TAB?(I0T),TAB^(I0T),TaB4(I0T)'TAB5(I0 
IT)' TAB6( IOT ) • CF( JOP) 
?043  FODMATI  1H0'4F10.5  ///) 

2044  FORMAT!////  6 X  '  4HSMAT • 6X ' 4HPR0D ■ 6X ' 4HACCP ' 4X ' 5HPF SAM ) 
">^45  FORMAT!/////  ?X'  6F10.4) 

AMT  s  SL  -  0.5»VPP 

CALL  TWOPI (AMT '  1 ) 

CALL  TEPDO(DEG,TAR3'^7'RAK'AMT) 

TAZ  s  RAK 

AMD  s  SL  -  PS 

CALL  TWOPI  (AMR  '  1  ) 

CALL  TERPO! DEG ' TAB5 ' 37 ' RAL ' AMR ) 

TA71  s  RAL 

ACS  s  RAK*CKKF  +  PAL 

CALL  TEpPOfDEG' TAB4' 37'SKO' AMT  ) 

TAZ2  s  SKO 

CALL  TEpPO<DFG'TAB6,^7'STOX,AMP) 

TA7^>  s  STOX 

APS  s  STOX*(1.0  +  (SKO*CKKFn 

PALT  s  SL  -  0»-5*VP 

CALL  TWOPI (FALT1 1 ) 

CALL  TERPOfDEG' TABl ' ?7'PKX'FALT ) 

TAZ4  s  PKX 

AC<  s  PKX*C5F 

CALL  TERP0(DEG'TAB2'^7'PKXA'EALT) 

TA75  s  PKXA 

APr  s  1.0  +  (C5F*PKXA) 

PRINT  2044 

PRINT  2043'  ACS'  APS1 ACk- 'APK 

PRINT  2045'TAZ'  TAZ1'  TAZ2*  TAZ3'  TAZ41  TaZ5 

RFTURN 

FNP1 

22 


PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SUBROUTINE  DAYXX( MONTH' DAY' ST T • TM ' DAYB' DAYM ' GRBS ' GRMS ' CP ) 


RPOO  FOPMATf  ////  5X'  8H  DAyB  s  'F5.0'10X'  8H  DAYM  s  'Fs.O  ///  5X'  8H  G 

1HBS  s  lF7.2'10X'  8H  GHMS  s  'F7.2  ///  5X'17H  ORIGINAL  T.M.  s  'F7.2' 

?10X'18H  rriRRECTFn  T.M,  s  'F7.2) 

GO  TO  (  207' 209' 21.1  '  213' 215' 21  7'219' 251  '  223' 225' 227*229)  'MONTH 

2^7  DAYP  s  DAY  207 

GO  TO  230  208 

209  DAYB  s  DAY  +  31.  ?09 

GO  TO  230  2  10 

21  1  OAYR  s  HAY  +  ^9.  211 

GO  In  2^0  2  12 

21^  DAYB  s  DAY  +  on.  213 

Go  Tn  230  214 

215  DAYB  s  DAY  +  120.  215 

GO  To  2^0  216 

217  DAYB  s  DAY  +  151.  217 

GO  To  2  30  218 

219  DAYB  s  DAY  +  181  .  219 

C,n    To  ?^n 

2=,i  HAYR  s  DAY  +  212.  221 

GO  TO  230  222 

22^  DAYP  s  DAY  +  243.  223 

GO  To  230  224 

225  DAYP  s  DAY  +  273.  225 

GO  TO  2^0  226 

227  DAYB  s  DAY  +  104.  227 

GO  TO  230  228 
2?9  DAYB  s  DAY  +  334. 

230  IF(STT.GT.12.0)  GO  TO  23  3  230 

rp  s  T^  +  ( ^TT#1 5 . 0 ) 

IF(CP.LE. 180.0)  GO  Tn  235 

?33  CP  s  TM  -  (24.oo  -  5TT)*15.0  233 

DAYB  s  DAYB  +  1.  234 

2i5  GPP?  s  CP/15.0  2^5 

GRMS  s  GPPS  +  12.00  236 

IF(ORMS.GT.24.0)  GO  to  240  237 

DAYM  s  DAYB  +  14.0 

GO  TO  242  239 
240  GRMS  s  GRMS  -  24.00 

DAYM,  s  DAYB  +  16.0 
242  PRINT  8800' DAYB 'DAYM' GRBS' GRMS' TM' CP 

PFTUPN 

FND 


SUPPOI  iT IMF  AZIM ( DX' VR ' A' VN' < ' J' OOMpy ' MFC > 

COMPUTE  MAJOR  OR  MINOR  AXIS  OR  COMPLJTF  NORTH  AND  EAST  COMPONENTS 
DIMFNSION  DX(J)'  VP(J)1  VN(J) 

do  in   i  s  1 ' J 
ion  VN  (  I )  s  (nx(i)  +  roMpv  -  a  1*0.0174^3'* 

Go  to  ( i ' 2  )  '  y 

1  DO  Aft  I  s  1  '  J 

56  VM ( I )  s  VR( I ) *COS ( VN ( 1 ) ) 
GO  TO  8  8 

2  I F ( NFC. EQ. 4. OR. NFC. Fo.i )  GO  Tn  109 
DO  7  7  I  si' J 

77  VM  (  I  )  s  VR(  I  )*SIN(  VN.  (  I  )  ) 

p  E  T  U  R  N 
88  DO  09  LM  s  1 ' J 

TR  s  VN(|_M) 

TA  s  VR(LM) 
09  VN ( LM)  s  SIGN ( TA • TR ) 
109  PFTUPN 

FND 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 


SUBROUTINE  ASTRO ( XYER ' DAYB • DAYM • GRBS ' GRMS ' JOBX ) 

DETERMINATION  OF  ORBITAL  FLEMENTS  FOR  OBSERVED  SERIES 

Common /COSTX/CXX  (30) ' OFX( 5) 
COMMON  /  LEAPYP  /  hay 'month 

olTO  FOPMATflHO'  3RHERPOP PAPAMFTFR  (I)  FXCFFDS  LIMITS  ) 

?]1]  FORMAT) 1H0' 15X'1HR' lOX'lHQ' 15X'9HU  OF  M(2)) 
2112  Fopmak  8X1  F10.3'  3X •  F10.3'  9X '  Fin. 4) 
7040  FOPMATI 1H1 ' 10F10.3////  IX • 10F10 .3/ ///  1X'10F10..3) 
TALL  ORBIT (XCFN'XSX'XPX'XHX'XPIX'XNX'OEX'T1 XYER '5) 
Do  30  NOE  s  1 ' 30 
TO  CXX(NOE)  s  0.0 

XCFT  s  XCEN  +  1.0 
DORY  s  0.25*(XYEP  -  XCFT ) 
AMI  s  XYFP  -  XCFM 
AMjT  s  AINT(DOBY) 
FA.pm  s  0.25*(AMI) 
FAPX  s  FARM  -  AINT(FAPM) 

I F ( FAPX.FO. 0.0. AND. MONTH. Fq.t )  GO  To  ^5 
1}      IF ( FARX. FQ.O.O. AND. DAYB. LF. 60.0 >  DAYB  s  DAYB  -  -j  .  o 
I  F  {  FAPX  .GT.O.  0  I  TiAYP  s  DAYP  -  1.0 

IF{ FARX. EQ. 0.0. AND. DAYM.LE. 60.0)  DAYM  s  DAYM  -  1.0 
IF ( FAPX. GT.O. 0)  DAYM  s  DAYM  -  l . o 
32  CXX(l)  s  XSX  +  129.38482032*AMI  +  13. ] 763Q6768* (DAY8  +  AMIT)  +  0.5 
lion] A^32*GRBS 

rXX(2)  s  XPX  +  40.66?46584*AMI  +  0 . 1 11 4040 1 6* ( DA YB  +  AMIT)  +  0.004 
1 641 R'4*GPBS 

CXX(3)  s  XHX  -  0.2387?498R*AMI  +  0 . 9 8 ^647' 2 P 8* ( DA Y R  +  AMIT)  +  0.04 
1  1  OA86^87*GPRS 

CXXI4)  s  XP1X  +  0.01717836*AMl  +  0 .000047064* ( DAYB  +  AMIT)  +  0.000 
1  ^01  961*GPRS 

CXX<5)  s  XPX  +  40.66246584*AMI  +  0 . 1 1 1 4040 16* ( DAYM  +  AMIT)  +  0.C04 
1641 834*GRMS 

CXX(6)  s  XNX  -  19.328185764*AMI  -  0 . 05295 39336* ( DAYM  +  AMIT)  -  0.0 
102?064139*GRMS 
IF ( JOPX.LF.15.AND. JOBX.GT.  0  )  GO  TO  40 
CALL  TWOPI  ( CXX '6 ) 
GO  To  4] 
T5  I F( DAY.FO.l . )  GO  TO  ?? 
GO  TO  31 

40  CXXI7)  s  XPX  +  40.66?46584*AM I  +  0 . 1 11 4040 1 6* ( DAYM  +  AMIT)  +  0.004 
i 641 B34*GRBS 

0XX(8)  s  XNX  -  19.3281  R5764*AMI  -  o  .OS295  I9^f>*  (  DAY"  +  AMIT)  -  0.0 
1022064139*GPBS 
CALL  TWOPKCXX'  8) 

41  AN  s   LXX(6')*0. 0174533 
AX  s   CXXI6) 

FYF  s    .9136949  -  0 . OT=  69 ?6*G0 S ( AN ) 
CXX(9)  s  ACOS(EYE)*57.?957795 

IF(rxX(9).LT.l7.0.OP.rxX(O).GT.T0.0)  PRINT  7110 
CIC-  s   CXX(  9)  *0. 0174533 
IFtCIG. FQ.O.O)  GO  TO  2^0 

I F ( AX.EO.O.O.OP.AX.Eo.l 80.0 )  GO  TO  ?30 
VXX  s   .0896831*SIN I  AN) /SIN(CIG) 
CXX  (]  )  s  ASIN(\/XX)*57. 2957795 

IFt AX. GT. 180.0. AMD.  r XX (  1  0  )  . GT  .  0 .0 )  CXX(IO)  s  -1.0*CXX(10) 
CVX  s  CXX  (  1  0)  *0.01  74"' 

FXX  s   .2067274*51 N< AN)* ( 1 .0  -  0.0194926*COS ( AN ) ) 
IF ( FXX.EO.0.0 )  GO  TO  202 

FZZ  s   .99798^2  +  0.20^7274*COS( AN )  -  0  .  00201 48*COS ( 2 . 0*AN ) 
IF  (  FZZ. FQ.O.O  ),  GO  JO  ?02 
EXFZ   s  EXX/EZZ 
IF(  FXFZ.r,T.^4^o.o)  Go  TO  ?0? 
CXX(ll)  s  ATA,\(  FXFZ  )*^7.2Qc;779c; 

I  F  (  AX.GT.  180.0. AND. CXX(  11  )  .GT.0.0  )  CXX'(ll)  s  -].0*CXX(11) 
rFy  s  CXX(  11  )  *0.0174^, 
?n?  vpx  s  SIN  (  CVX  )  /  (  roc  (  r\/x  )  +  o.3^4766/SIN(?.0*'"IG  )  ) 
rXX(12)  s  ATAN<VPX)**7. 2957795 

I F { AX.GT. ] 80.0. AND. CXX( 1 2 ) .GT.0.0)  GXXI12)  s  -1.0*CXX(12) 
pvo  s  CXX  (  12  )  #0.01  7i"' 

VPY  s  SIN ( 2-0*CVX) / (COS! ?.0*CVX >  +  0 . n7 26 ] 84 /S I N ( C I G ) **2 ) 
rXX(13)  s  ATAN(VPY)*57. 2957795 

I F ( AX.GT. 1 80.0. AND. CXXI 1 3 ) .GT.0.0 )  CXX<13)  5  -1.0*rxX(1T) 
PVCP  s  GXX ( 1 3 )*n.01 7AS^3 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

?m  dgx  s  rxx  (5)  -  rxx  (in 

tall  TWOP I ( PGX '  1 > 

xpn  s  PGX*n.0i7it;n 

rxx  (14)  s  PGX 

PAX  s  SIN(2«0*XPG)/(rOS(0.?*riGl*'*?/(6.0*SIN(0.**CTG)**2  >  ~  COS(2. 
T  0*XPG)  ) 

PXX  s  ATA.N(RAX)*57. 2957795 

UM?  s  2.0*(CXX(11)  -  CXX (10)) 

UL?  s   UM2  -  PXX 

IF(UL2.GT. 180.0)  UL2  s  UL2  -  180.0 

IF(UL2.LT. 180.0)  UL2  s  UL2  +  1R0.0 

CXX (15)  s  UL2 

7FC^     s  (  5.0*COS(CIG)  -  1.0)*SIM(XPG) 

?rr  s  (  7»0*COS  (CIG)  +  1  .0  )  *("OS  (  XPG  ) 

rA|_L  FITAN(ZFS'ZFr'OXX' 2 ) 

rPAV    s   0.5*(  MMi)   +  OXX   +  oo.n 

rALL  TWOPI  (  CRAV  ]  ) 

r  X  X  (  1  6  >  s  CRAV 

IF( J0BX.GT.15.AND. JORX.GT.O  )  Go  Tn  88 

PGXX  s   GXX ( 7  !  -  CXX ( 11  ) 

TALL  TWOPI  ( PGXX  '  ]  ) 

xx°r-  s  PGXX*0  .017453^" 

rxx(17)  s  PGXX 

PA,TX  s  rx.X  (  8  )  *n.ni  745^ 

FYFX  s  O.Pl^ft°4o  -  0.0^56Q26*COS ( BATX) 

CXX(?  )  s  ACnc  (  FYFX  )  *R7.?o^77Qc; 

U M 7 X  s  2.0*(CXX(11)  -  rxx(io)) 

fyit  s  rxx ( 20 ) *0.0] 74^?^ 

ZFXS  s  ( ^.0*COS ( FYI T )  -  1 . 0 ) *S I N ( XXPG ) 

ZFxc  s  <  7.0*rrt^  (  fy  it  )  +  i  .o  )  *cos(  xxpr-,  > 

CALL  FITAN<  ZFXc. '  ZFXr  '  oxxx  •  ?  ) 

CXXU9)  s   0.5*(UM2X)  +  oxxx  +  on.n 

rpaVX  s  C X X ( 1 9 ) 

CALL  TWOPI  I  <~PAVX  •  1  ) 

CXX(IP)  s  C  P.  A  V  X 
qg  no  33  i\]T  s  1  '  30 

nPr-p  s  CXX( NT )*7O000.  +  5.0 
3  3  CXX (NT)  s  AINT(nPOD)*n.000l 

dpjmT  ?04-O'(rxx(IAX)'IAX  s  I1  30) 

PRINT  2111 

PRINT  2112'  RXX'OXX'iim? 

RETURN 

END 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SUBROUTINE  ORBIT ( XCEN ' XSX • XPX • XHX • XP1X ' XNX ' OEX ' T » XYER » NNN ) 

DETERMINATION  OF  ORBITAL  FLEMENTS  FOR  BEGINNING  OF  CENTURY 

DIMFNSION  OFX  (NNN) 
g  fopmaT(]hO'  5F12.4'  5X'  F5.0'  2 F6 . 1 ) 

S  s  1^.1763968 

P  s  0.1114040 

XH  s   .985647? 

D]  s   .0000471 

XN  s    .0529539 

XCAM  s  XYER*0.01  +  0.001 

XCFN  s  AINT ( XCAN )*100.0 

T  s  -3.0 

YP  s  2.5 

GAT  s  1600.0 

no  io  jv  s  i' 30 

GP  s  (GAT)*n. 01/4.0  +  n. 00001 
FPX  s  APMGP) 
PP  s  AIMT(FPX) 
COL  s  FDx  -  PP 

if (col.lt.o.ot o:  go  ^  n 

IF (GAT. FO. XCFN)  GO  To  12 

YR  s  YR  -  1.0 

GO  Tn  9 
11  IF(  GAT. EO. XCFN)  GO  TO  1? 

9  GAT  s  GAT  +  100.0 
in  CONTINUE 
1?  T  s  (GAT  -  1900.0)*0.0l 

OFX(l)  s  270.4374222?  +  ^07.892*T  +  0  .  00 25 25*T**2  +  0 . OOOOO 1 89*T** 
13  +  YP*S 

OFX(2)  s  ^34.?2801944  +  1 09 . 0 3 220 5 5*T  -  0  .  0 1 n?4444*T**2  -  0.000012 
15*T**3  +  YP*P 

OFX(3)  S  279.69667778  +  0.76892^*1"  +  0  .  000  30  ?  5*T**7  +  YR*XH 

OFX(4)  s  281.2208333?  +  1.719175*T  +  0 . 00045 278*T**2  +  0.00000333* 
1T**3    YR*P] 

OEX(5)  s  259.1825333?  -  1 34 . 1 42 3972*T  +  0 . 00 2 1 0556*T**2  +  0.000002 
122*T**3  +  YR*XM 

CALL  TWOPI  (OFX  '5  ) 

XSX  s  OFX ( 1 ) 

XPX  s  OFX12) 

XHX  s  OEX(3) 

X^i  X  s  0EX(4) 

XNX  s  OFX(5) 

PRINT  8'  (OEX(JX)'JX  s  1'5)  '  GAT'  T'  YP 

RFTUPM 

END 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SURPOUTINF  TWOPI (  AUG'  10) 
C 

c 

C      ELIMINATE  NEGATIVE  Op  ANGLES  GPEATEP  THAN  360  DEGREES 
C 

DIMFNSION   AUG (JO) 
no  114  mo  s  1 ' 10 
1?0  IF ( AUG(MO)  ) H  ? ' 1  1 4'  1  1  1 
117  AUG(M^)  s  AUG(MO)+  160.0 

GO  TO  120 
111  IF (AUG (MO)  -  ^60.0  )  1 14' 114 ' 1 1 5 
116  AUG(MO)  s  AUG(MO)  -  160.0 

GO  TO  120 
114  COMTINUF 
RFTUPN 
END 

SUBROUTINE  FITANI  AUS ' AUC ' PT A ' JMAP ) 
C 
C 

C       PFTFRMINF  ARCTANGENT  AND  PLACE  IN  CORRECT  QUADRANT 
C 

RXG  s   AUS/AUC 

IF  DI VIDF  CH^rk-  14'  i  =, 

14  I F ( AUS) 16' 5 5 ' 17 

16  PTA  s  270.0 
GO  TO  88 

17  RTA  s  90.0 
GO   TO  88 

15  RTA  s  ATAN(BXG)*57. 2957795 
GO  Tn  ( 88  '  38 )  '  JMAP 

18  IF (AUS) 32 ' 32' 14 
12  IF ( AUG) 11 ' 55' 17 
n  RTA  s  180.0  +  RTA 

GO  TO  8  8 
-54  I F ( AUT ) 1R  '  1 7  •  R8 
15  PTA  s  RTA  +  180.0 

GO  TO  8  8 
M    RTA  s  PTA  +  360.0 

GO  TO  8  8 
R5  PTA  s   .0 
88  CONTINUE 

RFTUPN 

END 

SUBROUTINE  XMFAN ( VIN » FH' N ) 
C 

c 

C      MEAN  SERIES  FOR  NON-TIDAL  CURRENT  OR  MEAN  SEA  LEVFL 
C 

DIMENSION  VIN(N) 
15  FORMAT!////  1X'14H  SUM  OF  SEP  I ES • 3X • F 1 2 . 3  //  IX'  7HD I  V  I SOR ' 1 2 X '  I  1 0 
1  //  1X'14HMEAN  OF  SEP  1 ES •  5X • F  10  .  5  //  IX'IOHTOTAL  DATA'9X'I10  ) 

MA  s  N 

FH  s   .0 

DO  10  MA  s  I'M 

IF( VIN( NA) .EO.0.0)  MA  s  MA  -  1 
10  FH  s  EH  +  VIN(NA) 

BH  s  FH/FLOAT(MA) 

PRINT  15'  EH'  MA'  BH'  M 

RFTUPN 

FND 
h 
h 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SUBROUTINE  TERPO  (  ANG  «TXX  •  LI  D  •  ANSWER  •  AUGM  ) 

INTERPOLATE  TAPLFS  21 ' 22 ' 23 • 24 ' 25 ' AND  26 

DIMFNSION  ANG -(LIP)1  TxxtLI1^) 
1  FORMAT ( ///23H  JOB  TERMINATED 26HFAILFD  TO  INTERPOLATE  ANG. 

LIN  s  LID  -  1 

no  2  0    ION  s  l1 LIN 

I F ( AUGM. GT.ANGI ION) .AND. AUGM.LT. ANG( ION+1 ) )  GO  TO  ?8 

IF ( AUGM. FO. ANG ( ION > )  GO  Tn  ?g 

IF( AUGM.EQ.ANG( ION+] ) )  Go  to  27 

GO  to  20n 
26  ANSWFR  s  TXX ( ION ) 

GO  TO  201 
?7  AN.SWFR  s  TXX ( ION+] ) 

GO  Tr«  201 
?R  RFW  s  AUGM* 0.1 

WFP  s  REW  -  AINT(PFW) 

FPA  s  WFR*10.0 

PFT  s  APS ( TXX (ION) ) 

TFX  s  APS ( TXX ( ION+1 ) ) 

IF ( PFT.GT.TFX )  GO  TO  20 

IF(PFT.LT.TFX)  GO  TO  ^0 

IF(PFT.FO.TFX)  GO  TO  f\ 
79  DIFP  s  (PFT  -  TFX)*0.1 

IF ( TXX(  ION) .GE. 0.0. AND. TXX(  ION+i  )  .LT.O.P  )  DIFP  s  (PET  +  TEX)*0.1 

I F ( TXX(  ION )  .LT.n.O.AND.TXX (  ION  +  I  ) .GF.0.0  )  DIFP  s  (pET  +  TFX)*0.1 

DLFX  s  DIFR*FPA 

CHS  s  SIGN! DI FX  '  TXX  (  ION)  ) 

ANSWER  s  TXX(ION)  -  CHS 

GO  TO  201 
30  DIFP  s  (TFX  -  PFT)*0.1 

I  F  (  TXX  (  ION  )  .GE.O.O.AND.TXXI  ION+1  )  .LT.ri.f)  )  DIFP  s  (RET  +  TEX)*0.] 

I F ( TXX (  ION ) .LT .O.O.AND.TXX (  ION  +  l  )  .GE.0.0  )  DIFP  s  (PET  +  TEX)*0.1 

DIFX  s  DIFP#FPA 

CFK  s  SIGN ( DI FX' TXX ( ION+1 ) ) 

ANSWER  s  TXX(  ION)  +  CH^. 

GO  TO  2  01 
■>!      IF  (  TXX(  ION)  .GE. O.O.AND.TXX  (  ION  +  1  )  .GF.n. n  )  GO  TO  26 

I F ( TXX (  ION )  .LT.O.O.AND.TXX (  ION  +  1 ) .LT.0.0  )  GO  TO  26 

IF(TXX( ION) .GF. 0.0. AMD. TXX ( ION+1 ) .LT. 0.0)  GO  TO  32 

IF(TXX(  ION) .LT.O.O.AND.TXX (  ION  +  1  ) .GF. 0.0  )  GO  TO  3 2 

PRINT  1 

CALL  FXIT 
^2  DIFP  s  (TXXUON+1)  -  TXX  (  ION  )  )  *0.1 

DIFX  s  DIFR*FPA 

ANSWFP  s  TXX(ION)  +  niFX 

GO  TO  201 
?O0  CONTINUE 
20]  RETURN 

FN1^ 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SUBROUTINE  DETAZ (FIN ' AZO' NR' COMPV'VIN • XTEMP ' INC) 

DETERMINATION  OF  TRUF  AZIMUTH 
THIS  ROUTINE  IS  A  SPECIAL  PURPOSE  PROGRAM  FOP  USE  WITH  TIDAL  CURRENT 
ANALYSES  IN  THE  COA^T  AMD  GEODETIC  SURVEY 

DIMFNSION  FI N ( NR ) 'VIM (MR) ' XTFMP ( NP ) 
2  FOPMATf 1H1 ' 30H  AZIMUTH  OR  MEAN  FLOOD  DIP.  s  ' F7 ,2 ' 5X ' I  5 • 5X ' Fl 0 . 2 ) 
-3  FOP«AT(  ]H0  •  14HERR0R  IN  DFTAZ  ) 
ft  FOPMATflH]. '28HOPERATOR STOP  THIS  JOB) 

LMM  s 

SUMD  s  0.0 

AI7  s  AZO 

IF( AZO. EO. 360.0)  AIZ  s  0.0 

DO  10   KAZ  s  l'NR 

TDIP  s  FIN(KAZ)  +  COMPV 

IF( TDIP.LF.20.0)  TDIR  s  TDIP  +  ^60. 0 

IF(ABS(TDIP  -  AIZ  -  360.0). LE.  20.0)  GO  To  on 

IF(ABS(TDIR  -  AZO).LF.  20.0)  GO  TO  on 

IF(ABS(TDIR  -  AZn  +  ?60.0).LE.  20.0)  GO  To  on 

GO  In    too 
oO  LMM  s  LMM  +  1 

SUmD  s  SUMD  +  TDIR 
100  CONTINUE 

AZO  s  SUMD  /  FLOAT(LMN) 

IF  DIVIDF  CHECK  R'  5 
8  PRINT  6 
4  PRINT  3 

RFTUPM 
^  AZn  s  A I N T ( A Z O ) 

CALL  TWOPI  ( AZO' 1  ) 

PRINT  2'  AZO' LMN' SUMD 

IF ( INC.FO.00)  GO  TO  7 

CALL  TI LT (VIM ' XTFMP' MP) 
7  RETURN 

FMH 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

^UPPniiTINF  EDAT(MX' VIM'DIN'MA  «LX' ICOR 'CD'CV'LRD'XTFMP" INC) 

INFEP  DATA  CORRECTIONS  To  SERIES  FROM  MAGNETIC  TAPF 

DIMFNMON  ICOP(LX)  '  PIN (MA)  'VIN(MA)  • CDf  LX )  'CV(LX)  '  XTFMPf  MA ) 
COMMON/ FOP T/AO (  4  ) 'RO(i) 
<-,     format M  qX  '  A^  ) 

6  FOdmaT ( 60X ' F? .0 ' ?4X ' F4.2 ) 

7  foomAT ( 5X '  IS'^X'FS.O'SX'FR.^^X'I'i'^X'AM 

0  fodmaT ( ?4F3 .2  ) 
■\  0  FopMAT  (24F3.0  ) 

1 ?  FOPMAT (  1 1  OX ' F^ .0 ) 
M  FOPMAT IIQX'FR.O) 

REWIND  ° 

IF( INC.FO.OO)  GO  TO  14 

IF(LPD.EO.OO)  GO  To  t? 

RFAn(9" 51 ) ( XTFMPf IM) • IM  s  l'LPD> 
11  PFAD(9'P0)  (XTFMPf IP) ' IP  s  l'MA) 

OF'-1 1  MP  <=> 
14  IF  s  ] 

LAG  s 

LEN  s  1 

^O  QR  I  s  l'MX'LFN 

IF  (  [_PP>  QO  '  1  '  2 

2  IFf I    L"n) W] 

3  PFAD (P' p )   I  TRASH 
GO  TO  9  8 

1  LAG  s  LAG  +  1 

I F (  IF.GT.LX  )  GO  TO  1] 
LP  s  ICOPf IE) 
IFfLD.FO. I )G0  TO  4 
LDX  s  ICOR(IF)  -  1 
11   IFf  IF.GT.LX  )  I_PX  s  MX 
DO  40   IX  s  I ' LPX 
PFAD(P'AO)  DIN(LAG)'  VIN(LAG) 
I F ( IX.FO.LDX )  GO  TO  ^00 
LAO  s  LAG  +  1 
i^nc     oomTINUF 

I  F ( LDX.FO.MX )  GO  TO  07 
LFN  s  LD  -  I 
GO  TO  98 

4  CFAD(9*^)   ITRASH 
DIN (LAG)  s  CD(IF) 
VIM(LAG)s  CV( IF) 
XTFMP(LAG)  S  0.0 

PRINT  7'  I»  DIM(LAG)1  VlN(LAG)'  I  COP  (  I  F  )  '  I  TP  A  <^H 

LFN  s  1 

IF  s  IF  +  1 
08  CONTINUE 
0  7  IF(LAG.LT.MA)  GO  TO  q 

GO  TO  9  9 

8  LAO  5  LAG  +  1 

PFAD  01   (VIN(IN)   'IN  s  LAG'MA) 
PFAD  10'   (DIMflD)'  ID  s  LAG'MA) 
DO  3  0   MAG  s  LAG'MA 
^00  XTFMP(MAG)  s  0.0 
00  IF(INC.EO.l)  <"ALL  T  I  I..T  (  V  I  M  '  XTEMP  '  MA  ) 
IF ( INC.EO.l )  INC  s  0 
RETURN 
FMD 
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PROGRAM  FOR  HARMONIC  ANALYSIS  AT  TIDAL  FREQUENCIES-Continued 

SIIPPOIJTINF  TILT  (GIN'  TFM'L  IT) 

COPPECT  DATA  F01?  INSTRUMENTAL  TILT 

DIMFNSION  GIN  (LIT )  ' TFM ( L I T ) 
7    FOPMATf  1H0'  16  '  ^X  '?RHTIlT  fcXCEEDS  "*  5  DFGREF  LIMIT'3X'I3) 
no  10   IP  s  1 ■ L I t 
IF(TEM( IB) .EQ.0.0)  GO  TO  100 
MIX  s  IFI X ( TEM( IP ) ) 
NIX  s  MIX/5  +  1 
GO  TO  (100l10lI102'103,l'0  4,10  5'106,l07,110,110,110'lir)l110>  iNjX 

101  GIN (IB)  s  .03*GIN(IB)  +  GIN(IB) 

GO  TO  ICO 
10?  GIN (IP)  s  .06*6IN(IB)  +  GIM(IB) 

GO  TO  ]0O 
10?  GIN (IB)  s  .09*GIN(IB)  +  GIN (IB) 

GO  TO  100 
104  GIN (IB)  s  . 1?*GIN(IB)  +  GIN(IB) 

GO  TO  100 
1^5  g I N ( I P )  s  . 1 6  *G I N ( I B )  +  G I N ( I B ) 

GO  To  100 
106  GIN (IP)  s  . ?0*GIM(IB)  +  GIN(IB) 

GO  TO  100 
10  7  G I N ( I B )  s  . 2 5*G I N ( I B )  +  G I N ( I P ) 

GO  TO  100 
lOfl  GIN (IP)  s  . 51*GIN(IB)  +  GIN(IB) 

GO  TO  100 
10Q  GIN (IB)  s  .66*GIN(IB)  +  GIN(IB) 

GO  TO  mo 
110  PRINT  2  '  IB'MIX 
100  CONTlNUF 

RETURN 

END 

SUBROUTINE   T  FMP  (  XSFP  '  xTF^p  '  I  ) 

STOPF  SEPIES  FOP  FOUPIFP  ANALYSIS 
DIMENSION    XSFP(I)  •  XTFMP(I) 

no  i o  j  s  l '  i 

100  XTFMP(J)   s  XSFP (J) 
PFTUPN 
FNn 
SUB  POUT  INF  c-OPT  (  c>OPX  •  VX  '  KK  '  J  J  ) 

SORT  DATA  AND  FORM  NEW  SFPIFS 
DIMENSION   SOPXt  K<  )'  VX(  KK) 
COMMON /CONTT/NN 
1  DO  10   Is  l'KK'JJ 
NN  s  I/JJ  +  1 
]O0  VX(NN)  s  SORX ( I ) 
PFTUPN 

PNH 

SURPOUTIMF  ONFPK  AUX) 

TEST  FOR  180  DEGPEE  DIFFERENCE  BETWEEN  TWO  ANGLES 

IF(ABS(AUX)  -i  «0.0  )  1  70' 1  "70' if.T 

163  IF ( AUX )1 64 ' 1 64 ■ 1 65 

164  AUX     S     AUX     +     360. 0 
GO     TO     170 

IAS     Ally     s     AUX     -     360.0 
170    PFTUPN 
FNn 
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